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ABSTRACT 


CARRIER ESTIMATION USING CLASSIC SPECTRAL ESTIMATION 
TECHNIQUES FOR THE PROPOSED DEMAND ASSIGNMENT 
MULTIPLE ACCESS SERVICE 

BY 

BRADLEY JAMES SCAIFE, B.S. 

Master of Science in Electrical Engineering 
New Mexico State University 
Las Cruces, New Mexico, 1999 
Dr. Phillip L. De Leon, Chair 

In any satellite communication, the Doppler shift associated with the satel- 
lite’s position and velocity must be calculated in order to determine the carrier 
frequency. If the satellite state vector is unknown then some estimate must be 
formed of the Doppler-shifted carrier frequency. One elementary technique is to 
examine the signal spectrum and base the estimate on the dominant spectral com- 
ponent. If, however, the carrier is spread (as in most satellite communications) 
this technique may fail unless the chip rate-to-data rate ratio (processing gain) 


associated with the carrier is small. 



In this case, there may be enough spectral energy to allow peak detection against 
a noise background. 

In this thesis, we present a method to estimate the frequency (without knowl- 
edge of the Doppler shift) of a spread-spectrum carrier assuming a small processing 
gain and binary-phase shift keying (BPSK) modulation. Our method relies on an 
averaged discrete Fourier transform along with peak detection on spectral match 
filtered data. We provide theory and simulation results indicating the accuracy 
of this method. In addition, we will describe an all-digital hardware design based 
around a Motorola DSP56303 and high-speed A/D which implements this tech- 
nique in real-time. The hardware design is to be used in NMSU’s implementation 
of NASA’s demand assignment, multiple access (DAMA) service. 
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1 Overview 
1.1 Introduction 

In an effort to provide increased access to NASA’s Space Network (SN), a 
Demand Access Multiple Assignment communication scheme has been proposed. 
Under this scheme, users would have the option to communicate short information 
packets at low data rates on demand. This scheme is driven by the increased abil- 
ity of modern satellites (spacecraft) to detect error conditions on board the satel- 
lite [1]. Under current SN operations, communication services are pre-scheduled 
and the schedules often have significant delay and are not easily modified. The 
DAMA service is to be designed such that it operates independently of the current 
scheduled Multiple Access (MA) service. 

The SN consists of geostationary Tracking and Data Relay Satellites (TDRS) 
that operate as a virtual “radio frequency (RF) mirror” for data transmissions 
from/to low Earth orbiting (LEO) spacecraft communicating to/from ground sta- 
tions. These LEO spacecraft are in orbit about the earth and due to their relative 
motion to a receiving TDRS, a Doppler shift is induced in transmissions. Since 
the MA service is pre-scheduled, these Doppler shifts may be accounted for allow- 
ing the ground station to detect and receive the LEO signal. In the case of the 
DAMA system, where data transmissions are to be scheduled on demand, 
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the Doppler shift information may be unknown and thus a system must be de- 
signed to estimate (to within ±3 kHz) the carrier of the LEO signals. 

1.2 Comparison of Goddard and NMSU’s DAM A Proposals 

Two independent proposals have been offered to implement the DAMA ser- 
vice. The first has been proposed by N ASA/ Goddard Space Flight Center (GSF) 
and the second by New Mexico State University. The Goddard proposal intends 
to provide continuous tracking of all LEO satellites equipped with DAMA capabil- 
ity. Thus with the state vectors of all of these LEO satellites known, the Doppler 
shift of the carrier can be computed much like for the MA service. Whereas the 
Goddard proposal must maintain these state vectors so that the ground station 
receiver may demodulate, the NMSU proposal forgoes LEO satellite state vec- 
tor knowledge to simplify the required ground station hardware. In the NMSU 
proposal only a single element of the TDRS antenna array is used as a global bea- 
con [1], With a global beacon configuration, even a satellite that is experiencing 
alignment problems may transmit an emergency message to ground station users 
utilizing the DAMA communication system. As the state vector of the communi- 
cating spacecraft has been given up, the ground station will not know the position 
of transmitting LEO space vehicle and thus the Doppler shift cannot be accounted 
for [2]. Previous work has shown that the Doppler shift of a LEO spacecraft and 
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a TDRSS can vary by as much as ±50 kHz which is outside of the ground station 
receiver (GSR) tolerance of ±3 kHz [3]. 

The fundamental problem to NMSU’s proposal lies with estimating the carrier 
of the transmitted signal to within the tolerance of the GSR. Hardware must 
be developed that will provide a locking tone, accurate to within the ground 
station tolerance, in order for the ground station to demodulate. The problem is 
exacerbated by the nature of the DAMA signal. The DAMA carrier is to employ 
a Spread Spectrum (SS) scheme. SS signals tend to suppress spectral peaks of 
carriers and have responses that are spread out over a wider bandwidth. As we 
shall show below, this makes the proposed solution to carrier estimation more 
difficult. 

1.3 Proposed Solution 

The proposed solution to the carrier estimation problem described above is to 
employ classical Digital Signal Processing (DSP) spectral estimation techniques 
to estimate the carrier frequency. We shall employ a Discrete Fourier Transform 
(DFT) to generate magnitude squared spectral data from the received signal. The 
resolution of the DFT will be set to provide accuracy to within the ground station 
tolerance. A single iteration of this process will not be sufficient to protect the car- 
rier estimation from noise so we will average the magnitude squared data to limit 
the effects of noise. This process results in a periodogram. Having obtained the 
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periodogram of the received signal, we shall employ a frequency domain matched 
filter to maximize the Spectrum-to-Noise ratio (SPNR). The frequency domain 
matched filter will be predetermined based primarily upon the predefined SS fre- 
quency characteristics, such as processing gain (PG) and power, of the received 
signal. The results of the application of the matched filter, like those of the time 
domain equivalent, provide an optimal solution by enhancing the spectrum prior 
to searching. We then search the enhanced periodogram for a peak with which 
we base our carrier estimation. The accuracy of this technique will be shown 
to depend primarily on the PG of the DAMA carrier and the Signal-to-Noise 
ratio(SNR). 

1.4 Simulations, Test Data, and Theory 

At the core of the proof of concept for the proposed solution is a simulation 
model designed in Matlab. The simulation models the DAMA carrier against 
additive white Gaussian noise (AWGN), where we use AWGN to effectively model 
the white noise like spectrum of the MA service [2], Having approximated the 
TDRSS channel by its most critical feature(presence of the MA service) we add the 
DAMA carrier and simulate carrier estimation based upon the approach described 
in section 1.3. The simulation model provides estimation accuracy as a function of 
the SNR and the PG of the DAMA carrier. We will show that accurate estimation, 



to within the ground station tolerance, is achievable 80%-90% of the time for the 
given DAMA data rates and corresponding spreading rates. 

To verify the accuracy of the simulation model, we include the results of an 
experiment with test data gathered at NASA’s White Sands Complex (WSC). In 
this experiment, actual data vectors were streamed to a ground station transmitter 
for transmission to a TDRS (in orbit) and sent back to the ground station receiver. 
The parameters of the experiment were set such that we could observe several 
key issues with carrier estimation. The data gathered was processed with the 
proposed algorithm and compared to simulation results. The most significant 
conclusion of the experiment was that carrier estimation with the collected data 
was nearly equivalent to results obtained through simulation. The simulation is 
then recognized to accurately model the actual TDRSS channel. 

We have developed a theoretical analysis that leads to a rough approximation 
for carrier estimation accuracy. The analysis of carrier frequency estimation is 
based upon use of the DFT, and the description of carrier estimation accuracy as 
a random variable [4]. From this description and the use of various approxima- 
tions, we obtain an expression that describes the root mean square error (RMSE) 
between the actual carrier frequency and the estimate. The result is expressed as 
a function of SNR, data (chip) rate, and window type and length. Though the 
approximations break down in low SNR cases, in the higher SNR cases theory 
agrees with simulation results. 
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1.5 DAM A Hardware 


As described above, we seek to provide the ground station with an accurate 
carrier frequency based upon our estimation. We perform this estimation through 
the use of specifically designed hardware. The base of the hardware utilizes Mo- 
torola’s DSP56303EVM (EVM). The EVM utilizes Motorola’s DSP56300 core 
which is capable of 80 million instructions per second at 80 MHz and has enough 
available on-chip memory to implement the algorithm described above. As we 
will be required to sample the incoming signals at rates greater than that allowed 
by the EVM, we have integrated an 800 kHz 12 bit Burr-Brown ADS7810/19 
analog to digital converter (A/D) into the design. To interface the EVM and the 
A/D requires some additional logic and level translators that are implemented 
on an additional interface card. The card allows the EVM to control the A/D 
while allowing samples from the A/D to be passed directly into the memory of the 
EVM for processing. These three components, excluding some additional analog 
pre-processing and post-processing equipment, make up the core of the carrier 
estimation hardware. The hardware is designed to receive signals that have been 
filtered and frequency shifted to baseband, estimate the carrier, and then provide 
a locking tone to the GSR. The GSR will use this carrier estimate to demodulate 
the DAMA carrier. 


6 



The hardware has been tested with synthesized waveforms as well as actual wave- 
forms captured during the WSC experiment and performs as designed/required. 
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2 DAMA Project Description 
2.1 Current WSC Operations 


TDRSS was originally devised by NASA as an efficient means to control costs 
associated with providing a ground station for each satellite [5]. The concept of 
a space network was formed where users could transmit and receive all commu- 
nications through a common ground station. NASA operates TDRSS as a space 
network (SN) using it to provide customers with communication access to their 
Low Earth Orbiting (LEO) spacecraft. The SN consists of six geostationary Track- 
ing and Data Relay Satellites (TDRS) located 22,250 miles in orbit and a ground 
station located at the WSC (other operational ground stations exist as well) [5]. 
The function of a TDRS is to act as a virtual “RF mirror” through which commu- 
nication signals are relayed between user spacecraft and the ground station. An 
antenna array, located on each TDRS, is tuned by weighting antenna elements to 
provide a spot beacon to the spacecraft. This requires a unique weighting vector 
and associated signal processing equipment for each user spacecraft [1]. Two com- 
munication schemes are used by the SN to fulfill various communication needs [5] 
[3]: 

• multiple access (MA) at low data rates of 100 bps to 50 kbps operating in the 
S-band (2.1031 GHz - 2.1097 GHz forward service, 2.2845 GHz - 2.2905 GHz 
return) and using CDMA spread spectrum with a chip rate of 3 Mchips/s. 
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Figure 2.1: Illustration of Doppler Shift to DAM A Carrier 


• single access (SA) at up to 300 kbps operating in either the S-band or the K- 
band (2.0204 GHz - 2.1233 GHz forward service, 2.2 GHz - 2.3GHz return) 
using TDMA. 

The SN is able to provide 80% - 100% coverage for LEO spacecraft and is capable 
of simultaneously supporting 26 user spacecraft. 

TDRSS consists of geostationary satellites but the LEO spacecraft that use 
this system are not necessarily geostationary. We know that signals originating 
from a source moving relative to a TDRS will experience a Doppler shift [6]. It has 
been shown that the Doppler shift of these signals can be as much as ±50 kHz [3]. 
This is illustrated in Figure 2.1. The GSR normally maintains the state vector of 
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the satellite it is intended to communicate with and hence can simply calculate 
an estimate of the Doppler shift. Provided that the estimated Doppler shift is 
within ±3 kHz of the actual Doppler shift, the ground station can demodulate 
the received signal. In any system where one would forgo knowledge of the state 
vector of these satellites, the result would be that in general the GSR could not 
synchronize to the Doppler shifted carrier and thus the carrier need be estimated. 

The SN currently works under a scheduling process whereby a request for 
service must be made in advance (prescheduling) to utilize the SN. The scheduling 
delay often takes as much as 21 days for the request to be processed [3]. While the 
request can be serviced quicker in emergencies, the delay does not allow customers 
to react in near real-time to emergency situations that may arise with a user 
spacecraft. DAM A is a proposal that seeks to provide on demand communications 
between a user and their spacecraft without the need for prescheduling. 

2.2 NMSU’s Proposal 

The initial scope of the NMSU proposal is towards implementing a “911” 
service where satellites that have an gone into an error state may communicate 
this to the user when it detects such a condition. The eventual scope is to provide 
this service as a standard service to all DAMA capable spacecraft that require 
only low data rates with small data packets. Additionally it must be expanded 
to allow multiple access — or use by multiple users. For this thesis we assume a 
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single DAMA user at a time. The algorithm to be developed below is scalable to 
allow for multiple DAMA users at some point in the future. 

The NASA GSFC proposal seeks to implement the DAMA service by main- 
taining the state vector information for each user spacecraft that is DAMA capable 
by continuously tracking each of these spacecraft. As in the MA service, with the 
state vector of the spacecraft known, it is routine to estimate the Doppler-shifted 
carrier and provide this estimate to the GSR. In contrast, NMSU’s proposal gives 
up this state vector knowledge so that the DAMA ground station equipment is 
simplified. This leads to a problem with the Doppler estimation as it now must 
be estimated and supplied by means other than from the state vector. We pro- 
pose a solution to this problem with the algorithm to be developed below that 
will execute on the hardware that was also developed to provide the GSR with 
this Doppler estimate. DAMA is to be implemented with a SS BPSK modulated 
communication scheme like the MA service described above. However, there are 
certain restrictions that determine the parameters of the scheme. The DAMA car- 
rier is to be placed just inside the first upper TDRSS null as observed in Figure 
2.2. The carrier will be placed such that a maximal Doppler shift of +50 kHz will 
not place the carrier too near the null so that the rolloff of the TDRSS channel 
and other associated GSR equipment, which bandpass filters on the mainlobe, 
will not adversely affect carrier estimation. The signal in Figure 2.3 consists of 
the MA service and the DAMA signal and demonstrates the overall response of 
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Figure 2.2: DAMA Carrier Placement Against MA Spectrum 


the TDRSS system with the addition of the DAMA signal. We furthermore see 
the effects of sidelobe rejection of the TDRSS system. 

2.3 Spread Spectrum Fundamentals 

To discuss the operational parameters of the proposed DAMA carrier estima- 
tion, it will first be necessary to provide some fundamentals of SS communication 
schemes and definitions of important parameters. These parameters directly affect 
carrier estimation performance. 

We begin with a basic and widely used definition for SS systems: SS systems 
are distinguished by the characteristic that their signals consume a bandwidth 
greater than the information rate [7]. Though there are several different tech- 
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Figure 2.3: TDRSS MA Spectrum at IF from Actual Data 

niques for implementing the spreading, we focus on the technique known as spread 
spectrum by direct sequence (DS). DS spread systems implement a scheme where 
the information data is acted upon by pseudo-noise (PN) data, whose elements 
are referred to as chips, to produce a spread spectrum bandwidth (BW). The ratio 
of chips to bits is typically an integer and the chip rate is often much higher than 
the data rate. This ratio is defined as the processing gain (PG) where 

PG = A (2.1) 

Kb 

and R c is the chip rate in chips/s and Rb is the data rate in bits per second (bps). 
The PG also describes the ratio of chips/bit from which we see that each bit will 
be acted upon by PG chips through the use of modulo-2 addition. The PN code 
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Figure 2.4: 2047 PN Code Implementation 


sequence, Q, has the property that it approximates a white noise sequence and 
is periodic. C* is designed purposely such that 

CfCj = 8(i - j) (2.2) 

indicating orthogonality between PN codes of equal length but different “keys” . 
In practice, PN codes are only approximately orthogonal. The number of l’s and 
0 ! s, with C, G {0, 1}, differ by at most one. Many different techniques exist for 
the generation of these codes and we provide the 2047 PN code as an example. 
Though the 2047 PN code exhibits the qualities of white noise it is in fact periodic 
with period 2047. This PN code may be viewed as a primitive polynomial and 
implemented with a shift register as seen in Figure 2.4. 

In general the initial state of the shift register is a “key” and each key represents 
a different PN code that is orthogonal to other PN codes as described in (2.2). In 
this manner each PN code operates as an orthogonal basis function for each vector 
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of data. Multiple users are allowed in the same bandwidth precisely because each 
message is orthogonal to the other. 

For the MA service, R c = 3 Mchips/s. From theory it is known that the BW 
consumed by this modulation scheme follows the relation: 

BW a 2 R c (2.3) 

where BW is the bandwidth, and R c is the chip or spreading rate [8]. From (2.3) 
we observe that the MA service will occupy approximately 6 MHz of BW as is 
illustrated in Figure 2.2 and Figure 2.3. 

We will show that PG plays a large role in carrier estimation but first it is useful 
to see how PG will affect the BW of a signal. As we “spread” a carrier more and 
more (increase R c relative to R b ), the spectrum of the carrier will tend to spread 
out and flatten. This can be observed in Figure 2.5 below where we have spread a 
BPSK at various rates. Both PG’s in Figure 2.5 are relatively low however, it will 
be shown that low PG’s are required for accurate carrier estimation. Estimation 
of the carrier becomes more difficult at higher PGs since the carrier power is not 
concentrated over a small band of frequencies, which would result in a sharp 
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Figure 2.5: Spreading Effect of Processing Gains 


spectral peak, but rather distributed over a wider range. The peak in the SS case 
has less power and therefore noise may bury it (as intended for SS systems). 

2.4 Operating Parameters Description 

We now describe some of the communication parameters of the DAMA pro- 
posal. We are particularly interested in three parameters when dealing with digital 
communication systems: power, bit rate, and probability of bit error. For DAMA 
carrier estimation, we are not required to demodulate the signal and thus do not 
look at the probability of bit error. The first two, however, will affect our ability 
to to perform carrier estimation. 
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For the purposes of this thesis, we will describe DAMA power in two different 
ways. The first will be SNR in dB of the DAMA carrier-to-MA spectrum. We 
use this SNR definition when observing simulated work since the MA spectrum 
was modeled as AWGN. We can gain insight into this definition by observing 
Figure 2.3. The DAMA carrier will exhibit a peak against the passband of the 
MA service. Thus, we are interested in the DAMA power to that flat passband. 
The second way of describing power will be as the ratio C b /N a which is used in 
WSC operations. This is a measure of DAMA carrier power to the noise floor. 
Through observations of spectra (under typical conditions) collected at WCS, the 
mapping between the two power ratios was seen to be approximately 

C b /N 0 = 45dB » SNR = 2dB. (2.4) 

In this thesis, we shall use C b /N 0 to describe results of actual signals and SNR 
when describing simulation results. To avoid confusion when comparing the two, 
we shall map SNR to C b /N 0 . 

Since we intend DAMA to be a SS BPSK system, we must not only describe 
data rates, but also chip rates and therefore PGs. Due to the nature of the DAMA 
service, it has been proposed that R b = 1 kbps [3], The PG, as defined in (2.1), 
will be another parameter that is of primary concern since it determines the BW 
of the DAMA carrier and impacts the sampling rate. DAMA carrier estimation 
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accuracy will be shown to depend on PG to a large extent and therefore the chip 
rate will be a matter of investigation. 

Sampling rate is another important parameter. It was initially proposed that 
the PG would be set to PG = 100 [3]. With the data rate set as above, this 
implies a chip rate, R c = 100 kbps. This will exhibit a mainlobe width of 200 kHz 
by (2.3). From Figure 2.2 we can observe that we will need to account for another 
100 kHz due to Doppler shift. The total possible BW of the DAM A carrier is then 
300 kHz. We recall that to avoid aliasing while sampling the DAMA carrier, the 
DAMA signal must be bandlimited and be sampled at 

f s > 600kHz. (2.5) 

A commonly available, inexpensive 800 kHz A/D was found that matched the 
requirements and we thus chose f s = 800 kHz. This has implementation ramifi- 
cations that will be discussed below. 

2.5 Carrier Estimation Problem 

The problem with demodulating the DAMA carrier is the same as that for the 
MA carrier (though for the MA case the estimate is derived from the state vector 
of the satellite): the ground station is not capable of demodulating a signal if the 
error of the estimate is greater than ±3 kHz from the actual. As NMSU’s proposal 
will not keep track of state vector information for each user spacecraft, the carrier 
must be estimated reliably and efficiently and then passed to the ground station 
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receiver. NMSU’s DAMA proposal hinges upon accurate carrier estimation. We 
have developed an algorithm and hardware that performs this task and we shall 
describe the algorithm as well as the parameters for operation that will provide 
accurate carrier estimation. The hardware will provide a locking tone to the 
ground station receiver. 

2.6 Proposed Solution 


The proposed solution relies on classical spectral estimation theory with some 
modifications to improve performance. We begin by assuming that external analog 
hardware required to bandlimit and frequency shift to baseband the TDRSS and 
DAMA signals is available. Contained in the 400 kHz band (assuming f s ~ 800 
kHz) will be a portion of the TDRSS signal along with the entire DAMA signal. 

We next employ an averaged DFT (implemented with an FFT) which, if mag- 
squared values are computed, is also known as a periodogram. Since we assume 
AWGN with zero mean, the averaging has the effect of reducing estimation error 
variance of the carrier based on the DFT estimate of the DAMA spectrum. We 
express this as 




N - 1 

y: x[n + pN } e 

71=0 


j27rfcn 

N 


( 2 . 6 ) 


where P is the number of blocks in the average, N is the number of points in 
the block, n is the sample index, and k is the frequency index [9]. To obtain a 
frequency resolution that will enable us to estimate within the accuracy of the 
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GSR, we choose 



800 kHz 
3kHz 


« 267 


(2.7) 


but this would not allow us to use the radix-2 based FFT. We alternatively round 


to N — 512 for use with the FFT which yields 


f s mkHz 
N ~ 512 


1562.5 Hz 


( 2 . 8 ) 


This increases our increases our physical resolution beyond what is actually re- 
quired. 

From the result of the periodogram, we are left with an estimate of the spec- 
trum of the received signal. We estimate the carrier frequency by choosing the 
maxima of the periodogram. We assume that we are operating the DAMA service 
such that a peak will be observed in the average. We can improve the spectrum- 
to-noise ratio by utilizing a method from communication theory. It is known that 
the optimal solution for a receiver corrupted by AWGN is obtained by imple- 
menting a matched filter [7]. We will employ the frequency domain equivalent, 
spectral matched filter (SMF), which like its time domain analog is optimal and 
will maximize the SNR. In the time domain, a matched filter can be described by 
its impulse response 

h(t) = s(T - t) (2.9) 


with 


0 < t < T 
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where s(t) is the time reversed equivalent of the received signal. The matched 
filter is then convolved with the received signal yielding 


y(t) = [ s(t)s(T — t + r)dr. (2.10) 

Jo 

The same approach may be applied in the frequency domain where we have a 
SMF that is matched to the expected, frequency-reversed Power Spectral Density 
(PSD) of the DAMA carrier. Through a discrete convolution between the SMF 
and estimated DAMA spectrum, we may arrive at a desired optimal solution that 
maximizes the spectrum-to-noise ratio. The SMF is described as 

H(k) = X(N-k) , 0 < k < N (2.11) 

where X(N — k) is the frequency-reversed equivalent of the PSD of the DAMA 
carrier. The SMF is then convolved with the spectrum of the received signal as 

X sm f = X*H (2.12) 

where * indicates the convolution operator. The results of applying the SMF to 
the process will be described in detail in section three. We then form our carrier 
estimate with 

/ = ar8ma f U ' ) /, (2.13) 

where / is the estimated frequency of the carrier, N is the number of points of the 
FFT and f s is the sampling frequency. The discrete convolution in (2.12) smears 
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the number of points by 


L x smf = L H + M (2.14) 

where Lx am/ is the translated length of X sm f, L H is the length of the SMF (N = 
512), and M is the length of the periodogram (N = 512). We must subtract 
the maxima of the SMF, purposely located at N / 2 - or 256, from the index 
obtained by argmax(A r sm /). The SMF is generated in practice by generating a 
simulated DAMA carrier (without noise) with random data and averaging 1000 
periodograms based on the simulated signal. The simulated signal is set with the 
PG that we wish to test. A flowchart that describes the algorithm is shown in 
Figure 2.6. 
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3 Theory and Simulation 
3.1 Signal Description 

We have previously stated that the DAMA carrier will use SS BPSK and we 
begin with a description of the communication scheme. BPSK is a communication 
scheme where data is encoded into the phase of a sinusoid (carrier). The general 
equation for M-ary PSK is shown in (3.1) 

s(t) = Acos{cu c t + ip m {t)) (3.1) 

where A is the amplitude of the carrier, uj c is the carrier frequency and is the 
phase component. In general we may add more power to the sinusoid by increasing 
A but we assume A = 1 for the purposes of this thesis. ip m (t) may take on m 
distinct phases based upon the mapping of data to phase. For the binary case we 
use two distinct phases. Though we could choose any two distinct phases, in order 
to minimize error we typically choose the antipodal signaling as seen in Figure 
3.1 as it provides the optimal decision boundary [8]. With this configuration, 
the phase term xp(t) may take on values corresponding to ifr(t) 6 {0, n}. This is 
equivalent to modulating the sign of the sinusoid such that we may revise (3.1) to 

s(t ) = A d(t) cos(u; c f) (3.2) 

where d(t) € {—1, 1} and is dependent on the mapping of data. As explained in 
Section 2.3, the data is DS spread by a PN code. For purposes of the analysis that 
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Decision Boundary 



Figure 3.1: Signal Space Diagram for BPSK 

follows, we may look at SS BPSK as just BPSK at a higher data rate, namely the 
chip rate. We make this assumption since we are interested only in searching for 
a spectral peak of magnitude-squared data and not in demodulating the signal. 

3.2 Estimation Accuracy Theory 

The following arguments follow closely the work of Boaz Porat in his book, 
A Course in Digital Signal Processing. We have only extended his arguments to 
frequency estimation of a BPSK signal [4]. 

To begin, we make the following assumptions: s(n) is a BPSK signal as de- 
scribed in (3.2), S(e 3uin ) is the DTFT of s(n), w(n) is the rectangular windowing 
function, and v(n) is AWGN of zero mean. The received signal would then be 

x(n) = s(n)w(n) -t- v(n)w(n) 

taking the Fourier Transform of this yields 

N - 1 JV-1 

X( e ^ n ) = £ s(n)w(n)e~ 3un + £ v(n)w{n)e~ iun 

n=0 n=0 
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N-l N - 1 

= Y s(n)e~ Jun + Y v(n)w(n)e' 

n=0 n=0 

N-l 

= 5(e Ja ' n ) + Y v{n)w{n)e~^ n 

n=0 


\?um 


The PSD of s(n) is given by [8] 


Ps(e j n = E [sV**)] = ~ sine ( 


{u) — u; c )osf^ 

2tt t 


12 


(3.3) 


where we have transformed (3.3) from continuous time to discrete time, R is the 
data rate in bits/sec, and osf is the oversample factor in samples/bit. The maxima 
of the sine function as well as the maxima of the BPSK periodogram will occur 
at u = uj c hence 


N-l 


X(ei“ n )l =Uc = S(e> We ) + Y v(n)w(n)e>“' n 

71=0 

and taking the magnitude-squared of both sides, 


(3.4) 


j N-l \ ( N-l \* 

|X(u; c )| 2 = I 5(e , ' u ’ c ) + Y v{n)w{n)e~ iWc1x | f S(e JUJc ) + Y v{n)w{n)e~ iWcTl J 

V 71=0 / V 71=0 / 


N-l 


N-l 

= | S(e’ Uc ) + Y v(n)w(n)e -JWcn j ( S*(e ]Uc ) -f Y v(n)w(n)e 

V 71=0 / V 71 = 0 

= I 5 ( e ^)| 2 

N-l 

+23? 5(e jfu,c ) Y v(n)w{n)e-^ n \ 

n=0 

N-l N-l 

+ EE w(n)v(n)w(m)v(m)e 

n=0 m=0 


ju c n 


Taking the expectation of the magnitude-squared value yields, 
£[|A:(u. c )| 2 ] = E [|S(e7"')| 2 ] 
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+E 


23f? 


N - 1 

S{e? Uc ) £ v{n)w{n)e~i“' n 


n=0 


TN-\ N - 1 


El w(n)v(n)w(m)v(m)e Julci - n m ) 


Ln=0 m=0 

= P s (e ]UJn ) 


+25i 


iV— 1 


S(e jw ‘) £ E[v(n)}w{n)e~ j ^ n 


n = 0 


N-l W-l 

+ EE w(n)w(rn)E[v(n)v(m)]e JlAjc ^ n m 

n=0 m=0 


The first term is by definition as provided in (3.3) and is evaluated as 


= E 

The middle term evaluates to 0 since we have assumed E[v(n)] = 0. We further 
assume E[v(n)v(m)] — 7 v 8 {n - m), then 


E[\X(u> e 


| 2 l _ 
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Sj+ee w{n)w{m)^ v 5{n - m)e 2Wc{ - n m) 

n = 0 m=0 


A 2 


= i^ + ^E w2 ( n ) 


AR 
A 2 
4 R 


N - 1 

£ 

n=0 

+ 'y„N 


We next define an output SNR using (3.7) as the ratio between 


Ai 

SNR 0 = -^7 

7 iW 


(3.5) 


We define the input SNR as 

41 

SNRj = -2- (3.6) 

lv 
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where A 2 / 2 is the average power of a BPSK signal [8]. We delineate between 
output and input SNR to account for the application of the window function. 
Finally we define a window processing gain as 


W = SNR 0 
9 SNRj\N 

(3.7) 

— 

2 [ A 2 7„" 
N[ARj u Nf_ 

(3.8) 

— 

1 

RN 2 

(3.9) 


We next employ a rule of thumb which states 



NA 2 W 0 

(3.10) 


9 > 100 


lu 

- 

which is given in Porat as a requirement for the reliable detection and frequency 


estimation of a real sinusoid in noise [4]. We modify (3.10) for the BPSK 

case in 


terms of (3.5), (3.6), and (3.9) which yields 


— 

A 2 

(3.11) 


IV 
1— * 
o 
o 


In the development provided by Porat, the mean square error is then given as 

— 

F\ f f ] 2 ~ 24^0^ 

[Jc /cJ (2n) 2 A 2 D 3 

(3.12) 

_ 


where N a is the power density of the noise, J w is a window parameter [4], D is the 
measurement interval ( D = NT S , where T s is the sampling period), and A is the 
amplitude of the sinusoid. We further make the assumption that E[f c — f c ] = 0. 
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With this assumption and provided that (3.10) is true, the approximation in (3.12) 
is valid [4]. The approximation given in (3.12) can now be modified to yield 


E[f c ~ fc? 


24 N 0 J W 


(3.13) 


where we have applied (3.3) evaluated at uj = to c . The result of (3.13) is an 
approximation of the mean square error (MSE) of the BPSK carrier frequency 
estimate. We note that like the sinusoid in noise case, the approximation in 
(3.13) is only valid when we ensure that (3.11) is true. 

Equation (3.13) provides an approximation for the expected value of the MSE 
frequency estimate. We can use it to provide an approximation of the error in 
estimated frequency f c from a true frequency f c from using a DFT approach to 
estimation. It is recognized that with the many assumptions and the inclusion 
of a rule of thumb that (3.13) can only give a rough approximation to the actual 
MSE of the frequency estimation [4]. 

3.3 Comparison of Theory to Simulation 


Having obtained the MSE, we now show how (3.13) compares to simulation 
data. In the simulation we modeled a BPSK signal with a carrier frequency, / c , of 
178 kHz (typical of DAM A) and the sampling frequency, f s , set to 800 kHz (as in 
DAM A). The simulation performed an N — 512 point DFT on the BPSK signal 
and measured the MSE of the estimated frequency, / c , to the actual frequency f c . 
The results of the simulation are shown in Figure 3.2 where we have shown the 
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Figure 3.2: Comparison of Theoretical Curve Against Simulation 

simulation results versus the approximation developed in (3.13). The root mean 
square error (RMSE) is simply the square root of (3.13) and it is the expected 
value in Hz that f c is from f c as a function of SNR. This relates directly to the 
accuracy of estimating f c with f c to within the GSR tolerance and provides a 
means by which we can approximately determine the accuracy for a given SNR. 

We can conclude from the simulation that at high SNR’s (< 5 dB), the theo- 
retical curve and the simulation curve match, while at low SNR’s the theoretical 
curve becomes invalid (due to assumptions made during calculation) [4]. We 
further conclude that at high SNR’s (3.13) will provide a good estimate of the 
mean of the estimation accuracy and allow us to predict performance. At low 


30 



SNR’s the simulation curve approaches a condition where the noise is dominating 
the spectrum and the maximum peak is uniformly distributed over the 400 kHz 
band, yielding near random estimates. The theoretical curve is asymptotic to zero 
RMSE whereas the simulation curve will actually converge only to the difference 
between the true frequency and the selected DFT frequency — the difference is un- 
likely to be zero. This illustrates an important point regarding accuracy of carrier 
estimation: since we quantize to the DFT frequency points, there will most likely 
be an error irrespective of what SNR we are operating at. The maximum amount 
of this error assuming large SNR is given by 

max(/ c - /c) < ^ (3.14) 

This error can only be reduced by increasing N, the number of points of the DFT. 

3.4 Simulation Model Description 

Though we have developed a theoretical analysis culminating in (3.13), the 
approximation is not tight enough to provide proof of concept. We turn to a 
Matlab simulation model to provide a core proof of concept. We shall use this 
simulation model to explore the operational parameters of the algorithm. The 
simulation models the DAMA carrier against other MA users and AWGN. The 
simulation builds up a SS BPSK digital waveform based on the PG that we wish 
to test and then adds appropriate AWGN. It then estimates the doppler-shifted 
carrier frequency based upon the algorithm developed in section two. A record is 
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kept of carrier frequency estimations and those that fall within the tolerance of 
the ground station at WSC are counted as an accurate estimation. Likewise any 
estimation that is outside of the ground station tolerance is considered inaccurate. 
By performing the simulation 10,000 times for each SNR and desired PG, we 
obtain a plot that describes estimation accuracy. 

The simulation was originally written without the SMF and then rewritten to 
include the SMF process. The SMF improves estimation accuracy and allows for 
higher PG’s. The Matlab code for the simulation is given in Appendix A. 

3.5 Simulation Results 

We now provide results of simulation both in the SMF and non-SMF cases. 
The simulation results for the non-SMF are shown in Figure 3.3 and illustrate 
estimation accuracy as a ratio of accurate estimations to total estimations versus 
SNR. In Figure 3.3, we plot three curves for PG = 10, 20, and 100. The curves 
are generated through 10000 estimates per SNR and we perform 8 DFT blocks of 
N = 512 points. At an SNR = 2 dB, the results of Figure 3.3 demonstrate that 
estimation accuracy would be approximately 14%, 60%, and 88% for PG = 10, 20, 
and 100 respectively. NASA has placed the requirement that any implementation 
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Estimation Accuracy 



Figure 3.3: Family of Curves for Non-SMF Case 


should have as high a PG as possible. From Figure 3.3 it is clear that only the 
PG = 10 case is practical. 

Carrier frequency estimation may also be observed through the use of his- 
tograms which also relate information regarding the variance of the estimations. 
In Figure 3.4, Figure 3.5, and Figure 3.6 we show histograms for PG = 10, 20, 
and 100 respectively. Each is shown for SNR = 2. In each figure the actual carrier 
frequency is denoted by the center dashed line w T hile the GSR tolerance (±3 kHz) 
is shown by the outer dashed lines. 

The results with the SMF enhancement, as described in section two, are much 
better (Figure 3.7). The results shown in Figure 3.7 are generated as in the 
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Figure 3.4: Non-SMF Frequency Estimation for PG = 10, SNR = 2 dB 



Figure 3.5: Non-SMF Frequency Estimation for PG = 20, SNR = 2 dB 
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Figure 3.6: Non-SMF Frequency Estimation for PG = 100, SNR 



Figure 3.7: Family of Curves for SMF Case 
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non-SMF case with the exception that we have now included SMF processing. In 
Figure 3.7, we show estimation accuracy curves for PG = 10, 40, and 100. Here we 
see that estimation accuracy of the PG = 100 case has improved seven-fold. The 
case of PG = 40, a likely DAMA operating parameter, demonstrates accuracy 
at nearly 90% in a typical operating SNR region. As in the non-SMF case, we 
observe the histograms of Figure 3.8, Figure 3.9, and Figure 3.10 which are shown 
for PG = 10, 40, and 100 respectively. Each figure is shown at SNR = 2. From 
these figures it is clear that there is less variance and therefore more estimations 
within the GSR tolerance. 
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Figure 3.9: SMF Frequency Estimation for PG = 40, SNR = 2 dB 


2500 r 



Estimation Frequency (Hz) x 1 Q s 

Figure 3.10: SMF Frequency Estimation for PG = 100, SNR = 2 dB 


37 




3.6 Theoretical/Simulation Data Summary 


We have demonstrated results for a model that attempts to simulate DAMA 
carrier estimation in the presence of the TDRSS spectrum which is modeled with 
AWGN. We provide results indicating that accurate carrier estimation is possible 
with the algorithm developed in section two. Additionally, we have determined 
some operating points for the PG parameter. We have shown that we can ac- 
curately estimate the DAMA carrier 90% of the time with a PG = 40 at the 
typical SNR range. The simulation has provided a demonstration and we offer a 
loose theoretical approximation (3.13) to the estimation accuracy. In section four, 
we shall use data collected at NASA’s White Sands Complex (WSC) to further 
validate the simulation results. 

Appendix A includes all of the developed code for the simulation model in- 
cluding instructions for its use. Additional code is provided to perform data 
visualization. 


38 



4 WSC Data Collection Experiment 

4.1 Motivation of WSC Data Collection Experiment 

An experiment was devised to perform validation of the simulation model by 
collecting actual signals transmitted through TDRSS and collected at WSC. We 
subsequently processed them offline using the estimation algorithm. The funda- 
mental idea of this experiment is that we can observe the performance of the 
algorithm with actual DAMA signals. If the algorithm using actual signal data 
performs similar to the simulation under various PGs and at various SNRs, we 
may state with a degree of certainty that the simulation model is indeed accu- 
rate enough to predict carrier estimation. This also allows for the prediction of 
the performance of the algorithm if any of the operating parameters need to be 
changed. A side benefit is that we may also perform in-house testing using sim- 
ulated signals without the expense of interrupting critical TDRSS operations for 
testing. 

This section discusses the actual experiment performed at WSC, including 
details regarding setup, the processing of the data, and the conclusions drawn. We 
verify the simulation model’s results with actual TDRSS data. Additionally, the 
experiment points out the realistic operating boundaries which are an important 
part of the DAMA design. 
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Figure 4.1: Experiment Setup 


4.2 Data Collection Setup 

Two sets of equipment were used to collect test signals that were sent from 
the ground station and relayed off of a TDRS back to the ground station. The 
setup was as shown in Figure 4.1. The Transmit (Tx) equipment consisted of a 
computer equipped with a high speed PDMA32 Data Transfer Card (PDMA32). 
The purpose of the Tx equipment was to send a data vector to ground station 
equipment where it was BPSK modulated and transmitted to a TDRS. The data 
vector sent consisted of underlying data bits of {±1} with Rb = 1 kbps which 
was then spread by a spreading vector at a rate of R c = PG ■ Rb chips/s. Several 
of these data vectors, consisting of raw binary data and DS spread with a PN 
code (see section two), were generated in advance of the experiment and then 
used as a data source. The source code (wsands.m, wsands2.m, wsands3.m, and 
wsands4.m) is given in Appendix A. 
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The Receive (Rx) equipment consisted mainly of a LeCroy Digital Storage 
Oscilloscope (DSO) to capture IF signals and another computer to store captured 
signals. As the TDRSS channel is bandlimited to approximately 40 MHz, the 
DSO was set to sample at a rate of 100 MHz. This was the closest value to the 
Nyquist rate that the DSO was capable of sampling at. The DSO was capable 
of storing the captured data waveforms with either 50,000 samples or 100,000 
samples depending on the storage medium that was used and both sizes were 
collected. 

With the hardware setup described, we now turn out attention to the test set. 
The test set was established to test a variety of key parameters and determine the 
operational bounds of each. The key parameters are 

• Processing Gain 

• DAMA Carrier power-to-noise ratio ( Cb/N 0 ) 

• Placement of IF Carrier Frequency against TDRSS spectrum 

The PG, defined in (2.1), is the most critical of the three. Initial work with the PG 
indicated a PG of 100 could be used [3]. Subsequent simulations indicated that 
at this PG value the estimation accuracy was not reliable and lower values were 
investigated (see section three). It was included in the test set for completeness. 
PG’s of 10, 20, 40, and 100 chips per bit were the focus of the test set. 
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Figure 4.2: Placement of Test Frequencies Near TDRSS Null 

Three different carrier power-to-noise ratios were investigated. These ratios 
were 40, 45, 50 dB. These values were chosen as typical values based upon the 
experience of WSC staff [10]. It should be noted again that these values are 
described as a ratio of the carrier power to the noise floor in dB and not carrier 
power to MA spectrum as described in section three. 

The last test parameter is that of the carrier placement. The carrier frequen- 
cies were chosen such that a range of frequencies near the TDRSS null could be 
examined. Recall from section two that we intend to place the DAMA carrier at 
2.29 GHz which at IF is 2.9 MHz above the TDRSS center frequency. The range 
of test frequencies are observed in Figure 4.2. By examining a range of IF fre- 
quencies, we are allowed to effectively simulate a Doppler shift as well as examine 
the effects of the rolloff of the TDRSS channel. Comparing this to the TDRSS 
spectrum shown in Figure 2.2, we see that the upper range will be affected by 
the TDRSS rolloff. The full test set is shown in Table 1 . The full test set had to 
be reduced for logistical reasons to lessen impact on TDRSS operations, however, 
the reduced test set provides enough insight for simulation verification. 
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Table 1: Test Sets Captured 


C b /N a 

PG 

F c MHz 

40 dB 

45 dB 

50 dB 

10 

32.60 


* 


32.65 

* 

* 

* 

32.70 

¥ 

* 


32.80 


* 


20 

32.60 




32.65 

* 

* 

' * 

32.70 


* 


32.80 




B 

32.60 






* 


32.70 








* indicates data collected 


4.3 WSC Collected Data Processing 

We noted above that the captured waveforms were sampled at 100 MHz and 
stored as either 50,000 or 100,000 length vectors. In section two it was also noted 
that the bandwidth of the DAMA carrier plus maximum possible Doppler shift 
yielded a 300 kHz frequency search space and that the FFT resolution will be 
1562.5 Hz. Downsampling the captured waveforms to f s = 800 kHz was impracti- 
cal due to the problems of designing a narrow band decimation filter sharp enough 
for the rate conversion. Instead we seek a solution that will allow 7 us to perform 
a frequency search with a nearly equivalent FFT resolution. We see then that 

/ experiment = ~Tf — = 1525.8 Hz (4.1) 

N 65536 

provides a nearly equivalent FFT resolution to (2.8). In order to provide the 
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resolution in (4.1) and maximize the use of the short data sets, the vectors were 
either overlapped in the 100,000 sample case or zero-padded in the 50,000 length 
case. Without zero-padding the 50,000 length vector to 65536, the DFT does 
not meet the required resolution. The net result of this operation is a collection 
of 65536 length vectors with which to work with. These vectors were then used 
in a periodogram of eight data blocks chosen at random and processed with the 
simulation code that was used with a synthetic waveform (see Appendix A). This 
process was iterated 500 times, due to maximum data support, to form simulation 
bounds using WSC data. 

4.4 WSC Data Results Compared to Simulation 

The results of the processed data in comparison to simulation results are shown 
in Figure 4.3 below. The plot shown is for PG = 10 with the simulation repre- 
sented by the curve and the individual data points being the collected data across 
the various frequencies indicated. The plot indicates successful carrier estima- 
tion to within ±3 kHz as the ratio of accurate estimates to total estimates versus 
Cb/N 0 . We can make several observations from the plot. 

The first observation we can make is in regards to carrier power Cb/N 0 . We 
observe that at a Cb/N 0 = 40 dB, the estimation accuracy is much less than the 
simulation exhibits indicating a practical limit to C b /N 0 . As Cb/N a increases 
to 45 and 50 Cb/N a , estimates performed with actual WSC signal data match 
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Figure 4.3: Comparison of Simulated Results vs. WSC Captured Data 

those using synthetic data. This result indicates that we can estimate Doppler 
accurately provided that the carrier power is high enough for the range of carrier 
frequencies specified in Figure 4.2. This is an intuitive result as well since we 
would fully expect that a carrier with higher power will exhibit a more defined 
peak against the TDRSS spectrum. We also observe that the nominal carrier 
power provided by WSC will be sufficient as an operational bound. Previous 
work indicated that the DAMA carrier power and the MA carrier power can be 
relatively the same [3]. For the algorithm developed however we must have enough 
power so that the DAMA carrier exhibits a peak against the MA spectrum. Based 
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on collected data for the Cb/N 0 = 40 case, the DAM A carrier is lost in the MA 
spectrum and in general cannot meet reliability requirements. 

4.5 Conclusions of WSC Data Collection 

We have shown some important results with the WSC Data Collection analysis. 
The algorithm using actual TDRSS signals performs nearly the same as when 
using synthesized signals (see Figure 4.4). Therefore the synthesis of modeled 
signals is accurate to actual TDRSS signals. We rely heavily on the simulation 
to provide proof of concept and so the justification of the simulation is vital: the 
simulation is used to narrow in on the operating parameters that we should like 
to run at. Though we only have significant data for the PG = 10 case, we extend 
the results described above to higher PG’s with the assumption that the effects 
of spreading the signal even further in the frequency domain will cause neglible 
disparities between the simulation and actual operation. We can use simulated 
waveforms to implement additional, and important, real-time hardware tests. 

There remains one outstanding issue with regards to accurate carrier estima- 
tion. We have seen previously that the TDRSS channel begins to rolloff sharply 
near the location that we would like to place the DAMA carrier. Though the data 
did not support the investigation — at least in the spectral matched filter case — it 
is nevertheless true that if the DAMA carrier is nominally placed too close to the 
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Figure 4.4: Carrier Estimation Summary 


null, then the rolloff will adversely affect accuracy. We must also avoid placing 
the DAMA carrier in the MA passband. 

Overall the successful comparison of captured data to simulation data indicates 
that provided the DAMA carrier has sufficient power, we can explore additional 
operation bounds with a high degree of confidence. 


47 




5 Carrier Estimation Hardware and Software 
5.1 Motorola DSP56303EVM Description 

The DSP56303EVM (EVM) is the core hardware component of the carrier es- 
timation hardware developed. It executes the DSP specific real-time assembly 
language version of the carrier estimation algorithm (see Appendix B). The EVM 
is an evaluation module which is designed to be used in prototyping applications. 
As such, it offers a variety of interface options and configuration settings that 
make it adaptable to many different types of development projects. The board 
contains a DSP56303 24-bit digital signal processor that executes the assembly 
code routines. It also features a DSP56002 specifically for use in I/O functions 
through a JTAG/OnCE port. JTAG is a protocol that was developed to allow 
hardware and software developers to observe and manipulate hardware for trou- 
bleshooting purposes. The JTAG port for the EVM is primarily used with a PC 
and debugger software that allows a developer to load code into the EVM, single 
step through sections of code, and observe the contents of memory and registers. 
The onboard DSP56002 is not available for coprocessing code, but rather is used 
to control the I/O functions to the host PC. The EVM also contains 32K x 24-bit 
Fast Static Ram (SRAM) built with three banks of 32k x 8 bit SRAMS with 15 ns 
access times. An additional 64k x 8 bit Flash Programmable Erasable Read-Only 
Memory (PEROM) is provided for stand alone operation. The EVM also contains 
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Figure 5.1: DSP56300 Core System Block Diagram 


a 16 bit CODEC for sampling incoming analog waveforms and producing analog 
waveforms out from digital data [11]. 

The EVM is based upon the DSP56300 core which is shown in Figure 5.1. The 
EVM specifically uses the DSP56303 processor. The processor is capable of 80 
MIPS with an 80 MHz clock. It provides for backwards compatibility with 56k 
core code so that code written for earlier processors should function equally well 
on the newer cores. Due to a seven level instruction pipeline architecture, the 
DSP56303 is capable of effectively an instruction every clock cycle. It is based on 
the Harvard architecture so that it works with a dual memory structure. This 
architecture is particularly suited for parallel move instructions. This allows the 
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programmer to access two different memory structures at a time. For more on the 
DSP56303 please see [12]. 

The DSP56303 not only provides the processing power required to perform 
carrier estimation efficiently, but it is highly configurable and can interface to a 
large number of peripheral products. Though higher level languages (C) may be 
used to program the processor, it is most efficient when programmed in assembly. 
The assembly language that the DSP56303 uses is specialized to perform DSP 
tasks and includes features for that purpose. The Real-Time version of the carrier 
estimation algorithm is developed completely in assembly (see Appendix B). 

To create a program for the DSP56303, code is written in assembly language, 
and then assembled through the use of Motorola software tools to machine usable 
object code (.cld files). We may then use a software tool like Domain Technologies 
Debug-EVM to load the code into program memory of the DSP56303. With the 
Debugger software, we can utilize the JTAG port to single step through code and 
observe the results on calculations. 

For the DAMA project, as was described in section two, we are required to 
sample at a higher rate than the EVM can perform with the on board CODEC. 
We will use the on board CODEC in the generation of the locking tone to the 
ground station receiver. To sample the DAMA spectrum, we need to interface a 
high speed A/D. 
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5.2 Burr Brown 800 kHz A/D 


Earlier we noted we must search over a 300 kHz BW (200 kHz for the DAMA 
mainlobe, assuming no more than a PG = 100, and with ±50 kHz for maxi- 
mum Doppler shift) which was revised from an earlier estimate of 364 kHz. The 
Burr Brown 12 bit 800 kHz ADS7810/19 (BB A/D) was chosen since it has a 
sampling frequency of 800 kHz which gives a Nyquist interval of 400 kHz and 
because it is available in a convenient evaluation package for easy interfacing to 
the DSP56303EVM. The 12-bit samples of the BB A/D have a dynamic range of 
-72 dB which is sufficient for the TDRSS system. It was necessary to build an 
interface board to allow samples collected with the BB A/D to be passed to the 
303. 

5.3 A/D Interface Board 

The A/D Interface Board (ADIB) was developed as a Masters project by Tim 
Baggett [13]. Its purpose is to provide an interface to allow the EVM to com- 
municate with the BB A/D. The board was developed to allow' samples taken 
with the BB A/D to be passed into a peripheral (upper) memory location, which 
corresponds to the memory mapped I/O portion of memory. Though it is not a 
requirement, we may access the data samples with an efficient fast interrupt rou- 
tine which is preferable. It is also possible with the interface board to use direct 
memory transfers (DMA) though this was not implemented in the Real-Time ver- 


51 



sion of the carrier estimation algorithm. The board allows for user configuration 
of the specific memory location the samples will be written in to. The ADIB was 
further complicated by voltage level discrepancies between the BB A/D and the 
303. This was overcome through the use of zero wait-state level translators that 
adjust the output levels of the BB A/D to levels acceptable by the 303. 

5.4 Additional Hardware 

In an actual implementation of the DAMA carrier estimation hardware, there 
will be additional hardware required. Recall from the discussion of the DAMA 
project implementation, it was stated that WSC could provide the DAMA carrier 
at an IF of 32.65 MHz. The range of frequency values that a DAMA carrier can 
take, identified as 300 kHz, must be pre-filtered to bandlimit the signal containing 
the DAMA carrier from the MA waveform or aliasing will result. This is to be 
accomplished using an analog pre- filter set to a passband over the range if interest. 

With the signal now appropriately pre-filtered, the signal must be frequency 
shifted to baseband utilizing additional commonly available analog hardware. 
With this pre-processing accomplished, the DAMA carrier can be sampled and 
the digital signal processed to detect and estimate the carrier frequency utilizing 
the algorithm described in section two. Upon the algorithm’s determination of the 
DAMA carrier frequency, the CODEC on the EVM is utilized to provide a locking 
tone. By choosing a sample rate of 32 kHz for the synthesis of the locking tone, we 
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Figure 5.2: Locking Tone Generation 


need only frequency multiply by 25 to get the appropriate frequency value (< 400 
kHz). Afterward we frequency shift this tone back to IF for the GSR. This is illus- 
trated in Figure 5.2. The analog hardware required to implement the frequency 
multiplication and shifting is commonly available and an essential component in 
implementation. 

5.5 Software Description for Real-Time Carrier Estimation 

The carrier estimation code is designed to work with the EVM and the ADIB 
to sample the incoming signal, process the signal, and then produce a locking tone 
to the ground station. The core of the program is based upon the code supplied 
with the EVM named pass. asm. The purpose of the pass. asm is to initialize 
the EVM (CODEC and processor), then transfer samples to/from the CODEC 
from/to the processor. All signal processing is performed between transfers. This 
code is important because in most cases of programming a Motorola DSP56xEVM, 
the algorithm to be implemented will use pass. asm as the starting point. The 
carrier estimation code is slightly modified in that while it uses pass. asm as 
the core, it must set the EVM to properly receive samples from the BB A/D 
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and instead of reading sample values from the on-board CODEC, it will instead 
read sample values from an upper memory location as described above. Upon 
completing the initialization of the EVM (CODEC and processor) to work with 
the BB A/D, the carrier estimation algorithm begins. We begin by looking at the 
flowchart in Figure 2.6 which describes the algorithm. 

The initial step from the reset condition is to initialize memory and on-board 
CODEC and prepare the EVM for communication with the ADIB. The algorithm 
then begins by filling the sample buffer with values. Once it has achieved a 
full block of 512 samples, it applies a windowing function of the user’s choice. 
Typically either Rectangular or Hamming window coefficients are used. Included 
in this window function is an iteration scaling factor of 1/P, where P is the block 
number, which is nominally set to 0.125. This is included to scale for averaging. 
This has the added advantage of saving computation time since both windowing 
and average scaling are accomplished at once. 

From the flowchart we now see that the actual FFT is performed. The FFT 
routine is supplied by Motorola (see Appendix B) and produces a normally ordered 
Fourier transform, as opposed to a bit-reversed, complex result from a normally 
ordered input. The result is stored such that the real component of the result 
is in X memory and the imaginary in Y memory. We compute the magnitude- 
squared of the FFT and repeat this step averaging magnitude-squared data P 
= 8 times to form the periodogram. If we assume that x(n ) € 3? where x(n) 
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is the received, sampled signal, we need only search the first N /2 points due to 
Hermitian symmetry 

X(k) = X*(N - k), 0 < k < N - l. (5.1) 

This reduces the computational time and memory usage. 

Continuing with the flowchart in Figure 2.6, once the periodogram is com- 
puted, we convolve the result with the SMF and begin a search for the maxima 
of the enhanced spectrum. The result of this search is an index corresponding 
to the dominant spectral component. From this index a frequency may be found 
according to the relation 


/ = (5.2) 

where / is the estimated frequency [14], k is the index found from the search, 
L = N is the length of the SMF, N is the number of points of the FFT, and f s 
is the sampling frequency. A sine wave table is then used to provide a locking 
tone at 1/25 frequency that the ground station expects. As described above, this 
locking tone will be frequency multiplied up IF. The DSP56303 assembly code is 
given in Appendix B. 
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6 Conclusions and Future Work 


In this thesis we have described an algorithm to perform carrier estimation 
for the DAMA project. We have shown that the estimation accuracy depends 
primarily on the processing gain of the SS BPSK signal though some tertiary 
parameters also affect accuracy. We have developed a simulation model that 
provides proof of concept and provides DAMA operational parameters. Through 
the use of data collected at the WSC, we have verified the simulation model and 
established an operational point of PG = 40 with a corresponding accuracy of 
90%. 

In addition to the development above, we have developed Real-Time code 
based on the algorithm developed utilizing a Motorola DSP56303EVM. The entire 
project is implemented in hardware and functions equivalently to the simulation 
model. 

Future work might involve the formulation of a tighter theoretical bound as 
well the development of the analog hardware required for implementation. This 
work could be applied in other applications relating to communications between 
user spacecraft and commercial telecommunication satellites. In this scheme, ad- 
ditional carrier pre-shifting would be required. 
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A Matlab Simulation Code 
dama.m 

001 %******************************************************************* ****** 

002 X DAMA dama 

003 X 

004 '/, Description: This function allows simulation of the estimation of a DAMA 

005 X user’s Doppler-shifted carrier frequency. We assume a spread spectrum 

006 X system and the presence of other SS systems (modeled by AWGN) . 

007 7. 

008 X Programmer : Phillip De Leon 

009 7. 

010 X Creation Date: December 31, 1995 

011 7. Last Revision: Sept. 24, 1998 

012 7 . 

013 X Version History: 

014 7, 1.0 - frequency estimate of sinusoid in noise 

015 X 2.0 - frequency estimate of DAMA carrier in noise (spread spectrum) 

016 X 3.0 - better memory management, more efficient 

017 7. 4.0 - STFTs or average FFTs or MA filter smoothing of FFT bins 

018 X 5.0 - 800kHz sampling rate, 512-point FFT 

019 7. 6.0 - frequency estimate now over 612 blocks (of a longer signal). 

020 7. 7.0 - added frequency-domain matched filter to search 

021 X 

022 X Required subroutines: Use PLOT.DAMA.DATA.M to visualize data 

023 7 . 

024 7, Notes: 

025 X 1) For n MA users DAMA-to-MA power ratio is 

026 X SNR = (Rb_DAMA/Rc _DAMA)/(n*Rb_MA/Rc_MA) 

027 7. * ( 1000/ 10000) / (n* 10000/3000000) « 30/n. 

028 X For n = 5. SNR = 7.78dB; 

029 X 

030 X 2) This code is compute intensive!!! 

031 X 

032 X References: Digital and Analog Communication Systems by Couch p.359 

033 X 

034 7. Copyright (c) 1998 by Phillip De Leon, All Rights Reserved 

035 X* ************************************ ************************************ 

036 

037 clear: 

038 

039 7. 

040 X DAMA User's Parameters 

041 X 

042 Fc.DAMA = 164e3; X DAMA user’s carrier frequency (after shifting) in Hz 

043 main_lobe_BW = 200e3; X 200kHz main lobe BW centered on F.DAMA 

044 max_Doppler ■ 64e3; 7, maximum Doppler shift of +/-64kHz 

045 est.tolerance = 3e3; X must estimate carrier to within +/- est_tolerance 

046 

047 X 

048 X Communication System Parameters 

049 X 

050 Fs = 800000; X output signal sampling freq. (samples/s) 

051 Fc = Fc_DAMA + 0.6*max_Doppler; X assume DAMA carrier is +507, max 

052 X Doppler shifted (for testing purposes) 

053 

054 Rb.DAMA = 1000; X data rate for DAMA (bits/s) Fs/Rb_DAMA must be integer 

055 Rc.DAMA = 100000; 7. chip rate for DAMA (chips/s) 100K7 

056 samples.per.bit = Fs/Rb.DAMA; X must be integer 

057 samples_per_chip = Fs/Rc_DAMA; X must be integer 

058 Rb_MA = 10*Rb_DAMA; X data rate for MA (bits/s) 
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069 Rc.MA = 3000000; X chip rate for MA (chips/s) 3000K 

060 

061 X — 

062 X Frequency Estimation Parameters 

063 X 

064 X Estimation of Fc occurs over a bandwidth of... 

066 X Fe - max.Doppler - main_lobe_BW/2 < F.DAMA < Fc + max.Doppler + main_lobe_BW/2 

066 X 0 <* f <■ 328 

067 N « 512; X length data block to take FFT over 

068 window - 0; X if 0 assume rectangular window, 1 Hamming Window 

069 if (window) 

070 w = hamming(N) ; X Hamming window used in FFT 

071 end; 

072 

073 X— 

074 X Simulation Parameters 

076 X — 

076 number.of .estimates >30; X number of frequency estimates to perform 

077 X for each SNR typically 10000 

078 num.FFTs >8; X average 8 FFTs per estimate; 

079 snr « [-5]’; X vector of SNRs (in dB) to simulate typically from -2.2 to 4.7dB 

080 

081 X 

082 X Simulation 

083 X 

084 est.data * zeros(lengthCsnr) ,N/2) ; X malloc storage for estimation data 
086 


086 

087 

088 

089 

090 

091 

092 

093 

094 

095 

096 

097 

098 

099 

100 
101 
102 

103 

104 
106 
106 

107 

108 

109 

110 
111 
112 

113 

114 
116 
116 

117 

118 

119 

120 


for m = 1 : length (snr) ; X m is index for SNRs 


disp ( sprint f (’**» 
max.pos.est.error 
max.neg.est.error 
rand ('seed' ,0) ; 
randn( 'seed' ,0) ; 


Xd »**' , snr(m))) ; 
X reset variables 


SNR > 

- 0 ; 

= 0 ; 

X reset uniform generator seed 
X reset normal generator seed 


for n 

n; 


> 1 mumber.of.estimates; X n is index for estimate * 


sum.mag_X.sqrd 

5um.mag_S.sqrd 


zeros(N/2,l) ; 
zeros(N/2, 1) ; 


X malloc storage for accumulation of FFT data 
X malloc storage for accumulation of FFT data 


Synthesize SS signal 


Generate digital message 

data > round(rand(ceil(num_FFTs*N/samples_per_bit) , 1)) ; 


X generate data (bits) to modulate 


Digital BPSK Modulation 
start = 0; 

stop > samples.per.bit 
msg « [] ; 

for i « 1 : length (data) 
n > [start : stop] ; 
if ('data(i)) X if 
msg « [msg' 
else 

msg = [msg’ 
end; 

start • stop 


- l; 


0 multiply by -1 
(-l)*sin(2*pi»n*Fc/Fs)] 
X if 1 multiply by 1 
sin(2*pi*n*Fc/Fs)] ’ ; 


X build up message signal 
X build up massage signal 


+ 1 ; 


stop 

end; 


min((start + samples.per.bit - 1 ) ,num_FFTs*N) ; 


Direct Sequence Spreading 

PN.code » round(rand(ceil(num.FFTs*N/samples.per.chip) ,1)) ; 
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121 start = 1; 

122 stop = samples. per. chip; 

123 s = []; 

124 for i = 1 :length(PN_code) 

125 n = [start ; stop] ; 

126 if CPN.code(i)) X if 0 multiply by -1 

127 s ■ [s’ (-l)*msg(n) ’] ’ ; ‘/. build up spread signal 

128 else X if 1 multiply by 1 

129 s = [s’ msg(n)’]’; X build up spread signal 

130 end; 

131 start = stop + 1; 

132 stop = min( (start + samples. per. chip - 1) ,num_FFTs»N) ; 

133 end ; 

134 

135 X Compute FFT of BPSK-SS for Matched Filter Use 

136 S = fft{s(l:(0); 

137 mag.S.sqrd = real(S(l :N/2) .*conj (S(l :N/2) )) ; 

138 sum.mag.S.sqrd = sum.mag.S.sqrd + mag.S.sqrd; 

139 avg.mag.S.sqrd = sum.mag.S.sqrd/num.FFTs ; 

140 MF = flipud(avg.mag.S.sqrd) ; 

141 

142 for k «= lrnum.FFTs; X k is index on FFT number 

143 

144 X Add AWGN noise (MA users + channel noise) according to desired SNR to SS signal 

145 noise = randn(N.l) .* sqrt (cov(s)/( 10" (snr (m) /10) ) ) ; X ./ ? 

146 s.prime “ s(l+N*(k-l) :N*(k-l)+512) ; X Take length N buffer of samples from total. 

147 r = s.prime + noise; X carrier + noise 

148 r = r ./ sqrt(cov(r)) ; X scale to unit variance 

149 

150 X - 

151 X Window signal 

152 X 

153 if ("window) 

154 x = r; X no window => rectangular window 

155 else 

156 x = w . * r; X window with preset window 

157 end; 

158 

159 X = fft(x); 

160 mag.X.sqrd = real(X(l :N/2) . *conj (X(l :N/2) ) ) ; X due to symmetry need only scan lower half 

161 sum.mag.X.sqrd = sum.mag_X.sqrd + mag.X.sqrd; 

162 

163 end; X FFT loop 

164 

165 X 

166 X Frequency estimation 

167 X 

168 avg.mag_X.sqrd = sum.mag.X.sqrd/num.FFTs ; 

169 

170 X Frequency-Domain Matched Filtering 

171 MF.mag.X = conv(MF,avg_mag_X_sqrd) ; 

172 

173 if (0) 

174 k = f ind(avg.mag.X.sqrd == max(avg.mag.X.sqrd)) ; 

175 if (length(k) > 1) 

176 k - median(k) ; X choose the median k if there are several 

177 else 

178 k = k; 

179 end; 

180 end; 

181 if (1) 

182 k = f ind(MF_mag.X == max (MF.mag.X) ) ; 
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183 if (length(k) > 1) 

184 k • median(k) ; X choose the median k if there are several 

185 else 

186 k = k; 

187 end; 

188 k = k - N/2 - 1; X MF spreads range by 2 

189 end; 

190 

191 est_data(m,k) » est.data(m.k) +1; X tally estimate 

192 

193 end; X estimates loop 

194 save dama.data eat. data 

195 end; X snr loop 

196 

197 X 

198 X Save Data in MATLAB 4.x format 

199 X 

200 x = version; 

201 if (x(l)«“53) X running on MATLAB 6.x 

202 save -v4 dama.data Rc.DAMA Fs Fc snr est.tolerance N window number.of. estimates est.data 

203 else X running on MATLAB 4.x 

204 save dama.data Rc.DAMA Fs Fc snr est.tolerance N window number. of .estimates est.data avg.mag.X.sqrd 

205 end; 


60 



mf-bpsk_ss.m 


001 y,*************************»******«***«»»****»***********»**************»** 

002 */. MF.BPSK 

003 X 

004 */. Description: This function computes the matched filter for the PSD of the 

006 '/. BPSK-SS frequency-domain waveform. The matched filter coefficients are 

006 X stored in a file for use in the DAMA.M code. We may be able to use the 

007 7, theoretical PSD. 

008 */. 

009 X Programmer: Phillip De Leon 

010 X 

011 V. Creation Date: Sept. 25, 1998 

012 X Last Revision: Oct. 24, 1998 

013 */. 

014 X Version History: 

015 X 1-0 - MF for the BPSK-SS 

016 X 1.1 - New Vectorization Methods - will reduce run time. 

017 X 1.2 - Added Crash Recovery - ALL data is saved to temp. mat prior to 

018 */, next major algorithm step. 

019 ■/, 

020 X Required subroutines: 

021 X 

022 */. Notes: 

023 7. 1) N must be the same as in DAMA.M 

024 X 

025 Approximate Runtime: 

026 */. P200 vin98 : '22 minutes 

027 X Mac : Not performed 

028 y, Pxx Linux : Not performed 

029 X 

030 7, References: 

031 X 

032 X Copyright (c) 1998 by Phillip De Leon, All Rights Reserved 

033 7 t ************** *************************************************** ******** 

034 ON = 1; ‘/. Used to turn options on or off. 

035 OFF = 0; 

036 ERR = 0; X Set ERR to level achieved before crash - results 

037 X are saved in tamp. mat 

038 

039 if ERR == 0 disp( ’Operation Level One:’) 

040 X 

041 '/. Code Parameters 

042 X 

043 spreading. opt = ON; X ON to spread - OFF to just BPSK modulate. 

044 window « OFF; '/, if OFF assume rectangular window, ON Hamming Window 

045 

046 X 

047 7. DAMA User’s Parameters 

048 X 

049 Fc.DAMA “ 200e3; X DAMA user’s carrier frequency (after shifting) in Hz 

050 '/, Set to middle of DAMA frequency range. 

051 

052 y. 

053 X Communication System Parameters 

054 X 

055 Fs - 800000; X output signal sampling freq. (samples/s) 

056 Fc = Fc.DAMA; X assume unshifted DAMA carrier 

057 

058 Rb.DAMA = 1000; X data rate for DAMA (bits/s) Fs/Rb_DAMA must be integer 

059 Rc.DAMA = 50000; X chip rate for DAMA (chips/s) 100K? 

060 if rem(Fs .Rb.DAMA) '= 0 
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061 error( ’samples.per.bit MUST be an integer!!’) X Check to make sure that the 

062 X samples per bit is integer. 

063 end 

064 if rem(Fs,Rc_DAMA) *» 0 

065 error ( ’ samples.per. chip MUST be an integer!!’) X Check to make sure that the 

066 X samples per chip is integer. 

067 end 

068 samples.per.bit « Fs/Rb.DAMA; X must be integer 

069 samples_per_chip = Fs/Rc_DAMA; X must be integer 

070 

071 X 

072 X Frequency Estimation Parameters 

073 X 

074 X Estimation of Fc occurs over a bandwidth of... 

075 X Fc - mar_Doppler - main_lobe_BW/2 < F.DAMA < Fc + max_Doppler + main_lobe_BW/2 

076 X 0 <* f <» 328 

077 N = 512; X length data block to take FFT over 

078 

079 if (window) 

080 w « hamming (N ) ; X Hamming window used in FFT 

081 else 

082 w » ones(N,l); 

083 end 

084 

085 X 

086 X Averaging Parameters 

087 X - 

088 num.FFTs » 5000; X Average iterations - higher -> better resolution 

089 

090 save temp X End of level One - save work so far to temp. 

091 end 

092 

093 if ERR 1 I ERR «= 0 X Check Error Level. 

094 dispC’Operation Level Two:’) 

095 X 

096 X Generate digital message 

097 X 

098 M = num_FFTs*N*Rb_DAMA/Fs ; X number of bits to take PSD over 

099 

100 disp( ’Preallocating. .. ’) 

101 msg = zeros(M*samples.per_bit ,1) ; X Preallocate msg signal 

102 data = msg; X Preallocate data signal (-1.+1) 

103 s = msg; X Preallocate s - see variable description. 

104 dispCDone Preallocating...’) 

105 

106 msg * sin(2*pi*(0:M*samples_per_bit-l)*Fc/Fs) ; X BIG vector - not for the weak!! 

107 msg “ msg’ ; data * 

108 f ilter(ones(samples_per_bit , 1) ,l,upsamp((-l) . "(round(rand(M,l))) , samples.per .bit)) ; 

109 X Previous line generates a data vector with each 

110 X bit represented by samples_per_bit number of 

111 X samples. 

112 msg = msg.*data; X data vector is multiplied by the signal vector 

113 X to generate the message vector - BPSK modulated. 

114 

115 save temp X End of Level Two - save work so far to temp. 

116 end 

117 

118 if ERR «= 2 I ERR 1 I ERR *« 0; dispC’Operation Level Three:*) 

119 

120 X 

121 X Spread BPSK Signal 

122 X 
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123 if spreading.opt X If spreading option is on then spread signal. 

124 PN.code = f ilter(ones(samples_per_chip,l) ,1, ... 

125 upsamp(C-l) . ~(round(rand( ceil (num_FFTs*N/samples_ per .chip) , 1))) , . . . 

126 samples.per.chip)) ; 

127 s = msg.*PN_code; X s is the BPSK modulated spread spectrum signal. 

128 end 

129 

130 X Compute FFT of BPSK-SS for Matched Filter Use from averages . 

131 

132 X h = waitbar2(0, ’Generating Matched Filter Frequency Response...’); 

133 start = 1; 

134 stop = N; 

135 

136 save temp X End of Level Three - save work so far to temp . 

137 end 

138 

139 if ERR == 3 I ERR ■» 2 I ERR “ 1 I ERR ■» 0 

140 

141 disp( ’Operation Level Four;’) 

142 

143 sum_mag_S_sqrd = zeros(N/2,l) ; 

144 for 1 • l;num.FFTs X Generate averaged mag-squared 

145 X waitbar2(l/num_FFTs) ; X data from NN block FFTs. 

146 1 

147 n = [start : stop] ; 

148 temp » v.»s(n) ; 

149 S * fft(temp); 

150 mag.S.sqrd = real(S(l :N/2) . »conj CSC 1 : N/2) ) ) ; 

151 sum_mag_S_sqrd = sum_mag_S_sqrd + mag.S.sqrd; 

152 start = stop + 1; 

153 stop = stop + N; 

154 end 

155 

156 X close (h) 

157 avg_mag_S_sqrd = sum_mag_S.sqrd/num.FFTs ; X average. 

158 MF = flipud(avg_mag_S_sqrd) ; X "time reverse" data to prepare 

159 N1 = Fc*N/Fs; X for dama.mf code. 

160 X Record index of carrier frequency. 

161 

162 end 

163 

164 if window 

165 wtype = [’ha 1 ]; 

166 else 

167 wtype = [’re’]; 

168 end fname = [’pg’ int2str (Rc_D*MA/Rb_DAMA) wtype]; 

169 

170 mver = version if (mver(l) == ’5’) 

171 eval([’save -v4 ’ fname ’ MF N1 N ’ ] ) 

172 else 

173 eval([’save ’ fname ’ MF N1 N’]) 

174 end 

175 

176 

177 X 

178 X 

179 X Variable Description 

180 X 

181 X Name Description Purpose 

182 X 

183 X ON Program Control Used to turn ON/OFF Program Options 

184 X OFF Program Control Used to turn 0N/0FF Program Options 
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185 % 

ERR 

Program Control 

Used in disaster recovery 

186 X 

spreading. opt Program Control 

Used to turn ON/OFF spreading 

187 X 

window 

Program Control 

Used to turn QN/OFF windowing 

188 X 

Fc.DAMA 

Program Parameter 

DAMA Carrier Frequency - for the 

189 X 



purposes of this program it is used 

190 X 



to set the reference index. 

191 X 

Fs 

Program Parameter 

Sampling Frequency 

192 X 

Fc 

Program Parameter 

Carrier Frequency - in this program 

193 X 



it is set to the value of Fc.DAMA. 

194 X 

Rb.DAMA 

Program Parameter 

Bit Rate 

196 X 

RC.DAMA 

Program Parameter 

Chip Rate 

196 X 

samples.per. 

.bit 


197 X 


Calculated Parameter Each bit is represented by # of samples 

198 X 

samples.per. 

.chip 


199 X 


Calculated Parameter Each chip is represented by * of 

200 X 



samples . 

201 X 

N 

Program Parameter 

Length of FFT block 

202 % 

V 

Calculated Parameter window coefficients 

203 X 

num.FFTs 

Program Parameter 

Number of Iterations 

204 X 

M 

Calculated Parameter The number of bits required in the data 

205 X 



vector. 

206 X 

msg 

Calculated Vector 

The actual sinusoidal signal - later it is 

207 X 



BPSK modulated. 

208 X 

data 

Calculated Vector 

Binary data stored as (-1,1) 

209 X 

s 

Calculated Vector 

BPSK SS carrier 

210 X 

PN.code 

Calculated Vector 

The psuedo-random noise vector to cause 

211 X 



spreading. 

212 X 

start 

Program Parameter 

Used in indexing 

213 X 

stop 

Program Parameter 

Used in indexing 

214 X 

1 

Program Control 

Used in FOR loop 

215 X 

h 

Program Control 

Handle to figure 

216 X 

sum.mag.S.sqrd 


217 X 


Calculated Vector 

Holds sum of magnitude squared data from FFT. 

218 X 

S 

Calculated Vector 

Holds FFT data for current iteration 

219 X 

mag.S.sqrd 

Calculated Vector 

Holds magnitude squared data from FFT for 

220 X 



currant iteration. 

221 X 

avg.mag.S.sqrd 


222 X 


Calculated Vector 

Holds the average mag.S.sqrd vector. 

223 X 

MF 

Calculated Vector 

Frequency reversed - think of it as time 

224 X 



reversed data. 

225 X 

wtype 

Program Control 

Type of window employed for inclusion into 

226 X 



file name. 

227 X 

{name 

Program Control 

Filename 

228 X 

mver 

Program Control 

mver(l) holds matlab version. 

229 X 




230 X- 
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plot_dama_data.m 


001 %************************************************************************* 

002 X PL0T_DAMA.DAT A plot.dama.data 

003 V. 

004 % Description: This function provides a variety of tools for DAMA data 
00B X visualization using the file saved in DAMA.M 

006 X 

007 X Call Syntax: plot.dama.data 

008 X 

009 X Programmer: Phillip De Leon 

010 X 

011 7. Creation Date: December 31, 1995 

012 */, Last Revision: July 1, 1997 

013 7. 

014 X Required subroutines: 

015 X 

016 X Notes: Assume workspace (DAMA.DATA.HAT) has already been loaded. 

017 X 

018 X References: 

019 X 

020 'l Copyright (c) 1998 by Phillip De Leon, All Rights Reserved 

021 7>*********** ******************************************************** ****** 

022 

023 7clg; 

024 

025 selection = menu( ’Select a Plot’,... 


026 

’Estimation Ranges’ .... 

X 1) 

027 

’Estimation Histograms ’ , . 

.. 7 2) 

028 

’Estimation Accuracy’,... 

X 3) 

029 

’Exit’) ; 



030 

031 X 

032 7. Estimation range 

033 7. 

034 if (selection *« 1) 

035 i » f ind(est_data(l , : )) ; X find non-zero elements 

036 max_neg_est_error = ((min(i)-l) * Fs / N) - Fc; X 1 <= i <= N but need 0 <* i <= N-l 

037 max.pos.est. error = ((max(i)-l) * Fs / N) - Fc; 

038 plot([snr(l) snr(l)] , [max.neg.est.error max.pos. est. error] ) ; 

039 hold on 

040 for m = 2: length(snr) 

041 i = f ind(est.data(m, : )) ; X find non-zero elements 

042 max.neg.est.error = ((min(i) - 1) * Fs / N) - Fc; 

043 max.pos.est.error = ((max(i) - 1) * Fs / N) - Fc; 

044 plot([snr(m) snr(m)] , [max.neg.est.error max.pos. est. error] ) ; 

045 end ; 

046 plot([-100 100] , [est.tolerance est.tolerance],’:’) 

047 plot([-100 100] , [-est.tolerance -est.tolerance],’:’) 

048 hold off ; 

049 ylabelC ’Estimation error (Hz)’); 

050 xlabelCSNR (dB)’); 

051 if (window) 

052 t it le ([ ’Frequency Estimation Ranges (Rc=’ .sprintf ( ’7,d’ .Rc.DAMA) , ... 

053 ’, N=’ , sprintf ( ’%d’ , N),’, Hamming window)’]); 

054 else 

055 title( [’Frequency Estimation Ranges (Rc= ’ , sprintf ( ’7.d ’ .Rc.DAMA) , ... 

056 ’, N=’ .sprintf ( ’7d’ , N) , ' , Rectangular window)’]); 

057 end; 

058 [i.j] = f ind(est.data) ; 

059 Ymin = min(j); 

060 Ymax = max(j); 
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061 axis([(snr(l)-10) (snr(length(snr) )+10) l.l»min(Ymin,-est_tolerance) ... 

062 1 . l*max(Ymax,est.tolerance)] ) ; 

063 

064 % - 

065 X Estimation histogram 

066 X — 

067 elseif (salaction =■ 2) 

068 for m ■ 1 : langth(snr) 

069 bar( [0:N/2-l] . *(Fs/N) ,est_data(m, 1 :N/2)) ; 

070 hold on; 

071 plot ([0: M/2-1] ,*(Fs/N) ,est_data(m,l :N/2) , ’ + ’) ; 

072 plot([Fc Fc] , [0 1 . l*max(est.data(m,l :N/2))] , ’b: ’ ) ; 

073 plot([Fc+est_tolerance Fc+ast. tolerance] , [0 1 . l«max(est.data(m, 1 ;N/2) )] , ' r : ’ ) ; 

074 plot ( [Fc-est. tolerance Fc-est_tolerance] , [0 1 . l*max(est_data(m, 1 :N/2))] . ’r : ’ ) ; 

075 hold off; 

076 ylabel( ’Occurences ’ ) 

077 xlabel( ’Estimation Frequency (Hz)’) 

078 if (window) 

079 title( [’Frequency Estimation for SNR of ’ .sprintf (’’/.g’ , snr(m)), ... 

080 ’dB (Rc=* .sprintf ( ’Xd’ .Rc.DAMA) , N«’ , sprintf ( ’Xd ’ ,N) ,’ , Hamming window)’]); 

081 else 

082 X title ( [’Frequency Estimation for SNR of ’ .sprintf (’Kg’ , snr(m)), ... 

083 X ’dB (Rc*’ .sprintf ( ’Xd’ .Rc.DAMA) , N«’ .sprintf ( ’Xd’ ,N) ,’ , Rectangle window)’]); 

084 title([’Frequency Estimation for SNR of ’ .sprintf ( ’Xg’ , snr(m)), ... 

085 ’dB’]); 

086 end; 

087 axis( [0 . 96* (Fe-est_tolerance) 1 . 05« (Fc+est. tolerance) 0 1 . l*max(est.data(m, 1 :N/2) )] ) ; 

088 dispCHit any key to plot next histogram...’); 

089 pause ; 

090 end; 

091 

092 X 

093 X Statistical Accuracy 

094 X 

095 elseif (selection =» 3) 

096 lowar.index.bound » ceil(N»(Fc-est_tolerance)/Fs)+l X 1 <= lower_index_bound <= N 

097 X lower. indez.bound • 1; 

098 upper.index.bound « f loor(N*(Fc+est.tolerance)/Fs)+l X 1 <= upper. index.bound <= N 

099 X upper.index.bound « 211; 

100 

101 pause 

102 accuracy = zeros (length( snr) ,1) ; 

103 for m « 1 : length(snr) 

104 number.accurate_estimates°sum(est_data(m,lower_index_bound;upper.index_bound) ) ; 

105 accuracy(m) = number.accurate. estimates / number. of .estimates; 

106 end; 

107 plot(snr, accuracy, ’k’) 

108 xlabeK’SNR (dB)’); 

109 ylabel( ’Proportion of Estimates in Range’); 

110 title( ’Estimation Accuracy’); 

111 grid; 

112 else 

113 return; 

114 end; 

115 X plot.dama.data 
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wsands4.m 


001 y, White Sands Data Collection Code ver 4.0 

002 '/, Started on 6/26/98 

003 X Fill in comments . > Still need to comment 

004 •/, 

005 X Signal Vector Descriptions: 

006 X 1: Rb « 1000 bits/sec Rc = 100000 chips/sec 

007 ‘/. 2: Rb = 1000 bits/sec Rc = 20000 chips/sec 

008 X 3: Rb * 1000 bits/sec Rc = 10000 chips/sec 

009 X 4: Rb = 1000 bits/sec Rc • 40000 chips/sec 

010 X 

011 X Revision History: 

012 X 

013 X Ver 1.0: Baseline, one file generated per user input 

014 X Ver 2.0: Four "useful" files generated at once, no user 

015 ‘/, input. Some code cleaned up (vectorized). 

016 ’/, Ver 3.0: Hore vectorization, utilizes MATLAB’s filter 

017 X function. Vectors are pre-allocated . 

018 X Ver 3.01: Changed code so that # of bits may be individually 

019 */, selected. 

020 X Ver 4.0: To allov for larger data file generation, a menu 

021 ‘/, system is implemented so that only one data file 

022 ‘/, is generated. This compensates for memory 

023 ‘/. problems. Edit file for different spread/data 

024 ’/, rates . 

025 X Ver 4.01: "2047" NASA PN sequence took to long to generate 

026 '/, so the sequence is generated offline for faster 

027 X performance. This restricts resulting binary file 

028 '/, length to 2meg. 

029 X 

030 7. Last Update: 7/13/98 

031 ’/, Written by: Bradley James Scaife 

032 X* ******************* ********************************* ******************* 

033 timer ■ cputime 

034 preamblo.length * 5000; 

035 

036 choice = menu( ’Please Select File to be 

037 Generated : > , ’ lk.100 ' , ’ lk_20 ’ , * lk.10 ’ , ’ lk_40 ' ) ; 

038 

039 if choice == 1 

040 X Generate file one 

041 X,******************** 

042 

043 */, Message One Parameters 

044 X************************ 

045 '/msg.lenl = input (’Enter Message Length One (bits): ’); 

046 msg.lenl = 20e3; ’/.Restricted value - only change if spread vector code is altered. 

047 Rbl = 1000; Rcl = 100000; 

048 

049 ’/, Calculated Parameters: Message One 

050 '/,****************************»*»****** 

051 spread.f actorl = ceil(Rcl/Rbl) ; X Set as an integer. 

052 

053 ’/ Preallocate MSG1 vector 

054 X************************** 

055 msgl.init = zeros (msg.lenl , 1) ; msgl = 

056 zeros (msg_lenl*spread.f actorl , 1) ; 

057 Xspreading.vectorl = zeros (msg_lenl*spread_factorl, 1) ; 

058 

059 ’/, MSG 1 and Spreading Vector Generation 

060 X t **************************************** 
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061 rand( ’ seed’ ,0) ; 

062 

063 msgl. init ■ round(rand(msg_lenl,l)) ; 

064 

065 msgl = upsamp(msgl. init, spread. factorl) ’ ; 

066 msgl » filt«r(on«s (spread. factorl ,1) ,1 ,msgl) ; 

067 clear msgl. init 

068 

069 load svec 

070 

071 X Generate si 

072 X*************** 

073 disp( ’Generating si... 1 ) 

074 

075 si * 266*xor(msgl , spreading.vector) ’ ; 

076 

077 X Binary File Generation 

078 X************************* 

079 disp( ’Generating binary file...’) 

080 dispC Inserting preamble...’) 

081 disp( ’Inserting signal vector...’) 

082 

083 data.outl « [zeros (1 .preamble. length) si]; 

084 fid «fopen( ’ lk.100.bin’ , ’vb’ ) ; 

085 fvritadid, data.outl , ’int8’) ; 

086 fclose(fid); 

087 

088 elseif choice ■■ 2 

089 

090 X Generate File 2: Per definition in header above. 

091 X************* ************************ ************** 

092 Xmsg_len2 * input ( ’Enter Message Length Two (bits): ’); 

093 msg_len2 ■ 100e3; ^Restricted value - only change if spread vector code is altered. 

094 Rb2 * 1000; Rc2 » 20000; 

095 

096 X Calculated Parameters 

097 X************************ 

098 spread_factor2 * ceil(Rc2/Rb2) ; 

099 

100 X Preallocate MSG2 vectors 

101 X************************** 

102 msg2_init * zeros (msg_len2 , 1) ; 

103 msg2 > zeros(msg.len2*spread_factor2,l) ; 

104 Xspreading_vector2 «* zeros(msg_len2*spread_factor2,l) ; 

106 

106 X Massage and Spreading Vector Generation 

107 X****************************************** 

108 msg2_init = round (rand (msg.len2, 1)) ; 

109 msg2 » upsamp(msg2_init .spread.f actor2) ’ ; 

110 msg2 » filter(ones(spread.factor2, 1) ,1 ,msg2) ; 

111 clear msg2_init 

112 

113 load svec 

114 

115 X Generate s2 

116 X*************** 

117 disp( ’Generating s2...’) 

118 

119 s2 « 255*xor(msg2, spreading.vector) ’ ; 

120 

121 X Binary File Generation 

122 X********** *************** 
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123 disp( ’Generating binary file,..’) 

124 disp( ’Inserting preamble...’) 

125 disp( ’Inserting signal vector...') 

126 

127 data_out2 = [zeros(l .preamble. length) s2] ; 

128 

129 fid = fopenO lk_20.bin’, ’Mb’); 

130 fvrite(f id,data_out2 , ’ int8 ’ ) ; 

131 fclose(fid) ; 

132 

133 

134 elseif choice »= 3 

135 X Generate File 3: Fer definition in header above. 

136 x********************************** ****** *********** 

137 Xmsg_len3 = input( ’Enter Message Length Three (bits): ’); 

138 msg_len3 = 200e3; ^Restricted value - only change if spread vector code is altered. 

139 Rb3 ■ 1000; Rc3 = 10000; 

140 

141 X Calculated Parameters 

142 X************************ 

143 spread.f actor3 = ceil(Rc3/Rb3) ; 

144 

145 X Preallocate MSG3 vectors 

146 X****** ************** ****** 

147 msg3_init = zeros (msg_len3 , 1) ; 

148 msg3 = zeros (msg_len3*spread_factor3,l) ; 

149 spreading_vector3 * zeros (msg_len3*spread_f actor3 , 1) ; 

150 

151 X Message and Spreading Vector Generation 

152 X****************************************** 

163 msg3_init « round(rand(msg_len3 , 1) ) ; 

154 msg3 *upsamp(msg3_init,spread_factor3) ’ ; 

155 msg3 *> filter (ones (spread.f actor3 , 1) ,1 ,msg3) ; 

156 clear msg3_init 

157 

158 load svec 

159 

160 X Generate s3 

161 X* ************** 

162 disp( ’Generating s3...’) 

163 

164 s3 * 265*zor(msg3,spreading_vector) ’ ; 

165 

166 X Binary File Generation 

167 X************************* 

168 disp( ’Generating binary file...’) 

169 disp( ’ Inserting preamble...’) 

170 disp( ’ Inserting signal vector...’) 

171 

172 data_out3 = [zeros (1 .preamble.length) s3) ; 

173 

174 fid = fopenC ’ lh_10.bin’ , ’wb’ ) ; 

175 f vrite(f id,data_out3 , ’ int8’ ) ; 

176 fclose(fid); 

177 

178 else 

179 X Generate File 4: Per definition in header above. 

180 X*************************************************** 

181 Xmsg_len4 = input (’Enter Message Length Four (bits): ’); 

182 msg_len4 * 50e3; XRestricted value - only change if spread vector code is altered. 

183 Rb4 = 1000; Rc4 = 40000; 

184 
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185 X Calculated Parameters 

186 %*••*•*•*••••••*****••*** 

187 spread_factor4 « ceil(Rc4/Rb4) ; 

188 

189 X Preallocate MSG4 vectors 

190 %*••*****•••*•***••*••****• 

191 msg4.init « zeros(msg_len4,l) ; 

192 

193 msg4 - zeros(msg_len4*spread_factor4, 1) ; 

194 % spreading. vector4 « zeros (msg_len4*spread_factor4 , 1) ; 
198 

196 X Message and Spreading Vector Generation 

197 %«************************************»«*** 

198 msg4_init « round(rand(msg.len4,l)) ; 

199 

200 msg4 « upsamp(msg4_init,spread_factor4) ’ ; 

201 

202 msg4 = filter (ones (spread_factor4 , 1) , 1 ,msg4) ; 

203 clear msg4.init 

204 

205 load svec 

206 

207 % Generate s4 

208 %*************** 

209 disp( 'Generating s4...’) 

210 

211 s4 « 26S*xor(msg4 .spreading.vector) ’ ; 

212 

213 y. Binary File Generation 

214 %»*******»**»******»***»** 

215 disp (’Generating binary file...') 

216 dispC Inserting preamble...’) 

217 disp( ’Inserting signal vector...’) 

218 

219 data_out4 = Czeros(l .preamble. length) s4] ; 

220 

221 fid - fopent ’ lk_40.bin’ , ’ wb’ ) ; 

222 fvrite(f id.data.out4, ’int8’) ; 

223 fclose(fid); 

224 

225 end 

226 

227 dispCAll Done :) ’) 

228 run.time = cputime - timer 
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cap_mf.m 


001 7 Title: 

002 7 cap.mf.m 

003 7 

004 '/, Purpose: 

005 7 The purpose o f this code is to generate estimation data similar to 

006 % dama.mf with the exception that this code is tailored to use data 

007 7, captured during the White Sands Complex (WSC) experiment. 

008 % 

009 7, Revision: 

010 7. 1.0 — 11/20/98 

011 X 

012 % Revision History: 

013 7. none => baseline 

014 X 

015 7 Author: 

016 7. Brad Scaife 

017 I. 

018 7. Notes: 

019 7 This code requires as input data that has been prepared by wsco.d2.rn which 

020 7 prepares the captured data and generates a matrice of usable vectors (see 

021 7 vsco_d2 comments for more information). This code also relies on an index 

022 7 matice that is generated by in.prp.m. In both cases these files have already been 

023 7 prepared and may simply be loaded following prescribed naming conventions. In 

024 7. the event that these data files have been lost though it was thought helpful to 

025 7 comment on how to regenerate them. 

026 7 

027 7 Average runtime: 

028 7 "8 minutes 

029 7 

030 7******************************************************************************************** 

031 clear all; 

032 

033 ON « 1; OFF = 0; 

034 

035 7 Code Options 

036 7::::::::::::::: 

037 window * OFF ; 7 Use to toggle use of hamming window 

038 matched.f ilter = ON; 7 Use to toggle use of matched filter 

039 insight = OFF; 7 Use to gain instantaneous frequency estimation 

040 

041 7 Code Parameters 

042 7 ::::::::::::::::: : 

043 Fs = 100e6; 7 DO NOT CHANGE - this was the sample rate used. 

044 7snr =1.5 7 As per dama.mf definition - not wsc definition. 

045 est .tolerance = 3e3;7 +/- range of acceptable error in Hz. 

046 N = 65536; 7 FFT block length. 

047 number.of .estimates = 1000; 7 * of estimation attempts. Don’t modify!!! 

048 fft.avg =8; 7 Number of FFT’s to average over - changing this will 

049 7 require changes to code mentioned in documentation above 

050 7 as well as the rer unn ing of this code. 

051 

052 alpha =2; 7 Constant used to obtain consistent data points - does not 

053 7 add to statistical meaning nor detract from it. 

054 f.lo = 32.50e6; 

055 f.hi = 32.90e6; 

056 

057 lower .bound. index = ceil(f _lo*N/Fs)+l ; 

058 

059 start ing.bound = lower.bound.index; 

060 
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061 upper.bound. index « f loor(f_hi*N/Fs)+l ; 

062 Rb.DAMA - le3; 

063 

064 X Program Flow 

066 X: ::::::::::::: : 

066 if window 


067 

wtype « 

[>h*]; 



068 

else 





069 

wtype « 

t*r'] 




070 

and 





071 






072 

manu.sel « 

menu( 

Select Data 

to Process ; ' , . 

. 

073 

’SF*10, 

SNR1, 

F2 ’ , . . . 

X Selection 

1 

074 

’SF«10, 

SNR1, 

F3 * , . . . 

X Selection 

2 

075 

’SF«10. 

SNR2, 

FI' ... . 

X Selection 

3 

076 

•SF=10, 

SHR2, 

F2' , . . . 

X Selection 

4 

077 

’SF»10. 

SNR2, 

F3’ , . . . 

X Selection 

5 

078 

•SF»10, 

SNR2, 

F5\... 

X Selection 

6 

079 

■SF»10. 

SNR3 , 

F2 ’ , . . . 

X Selection 

7 

080 

’SF=20, 

SNR2, 

F3’) 

X Selection 

8 

081 






082 

if manu.sel 1 





083 data.f ilename « [ , c:\researeh\code\matlab\whiteB*2\uscmafl\snrl_f2' wtype] ; 

084 aval([’load 1 data.f ilaname] ) ; 

086 num.col ■ size(x,2); 

086 index.f ilename » [’e:\research\eode\matlab\whites"2\cmb8_’ int2str(num_col)] ; 

087 evalCC’load 1 index.f ilanama] ) ; 

088 itarations » choose (num.col , 8) ; 

089 Xsnr « -6; 

090 snr « 40; 

091 Rc.DAMA » 10a3; 

092 Fc - 32 .65a6 ; 

093 

094 alsaif menu, sal •= 2 

096 data.f ilanama = [’ c : \research\code\matlab\whites'2\wscmat'l\snrl_f 3 ’ wtype]; 

096 avalCC’load ’ data.f ilanama] ) ; 

097 num.col » siza(x,2); 

098 index.f ilanama « ['c:\research\code\matlab\whites-2\cmb8_' int2str(num_col)] ; 

099 eval(['load ' index.f ilename] ) ; 

100 itarations = choose (num.col ,8) ; 

101 Xsnr = -5; 

102 snr = 40; 

103 Rc.DAMA - 10a3; 

104 Fc « 32.70e6; 

105 

106 

107 alsaif manu.sel »= 3 

108 data.f ilanama * [’c : \research\coda\matlab\whites"2\wscmat"l\snr2_f 1 ’ wtype]; 

109 aval([’load ’ data.f ilename] ) ; 

110 num.col » size(x,2); 

111 index.f ilaname » [ , c:\research\code\matlab\whites"2\cmb8.’ int2str (num.col)] ; 

112 eval([’load ’ index.f ilename] ) ; 

113 itarations • choose(num_col,8) ; 

114 Xsnr « 1.6; 

115 snr * 45; 

116 Rc.DAMA - 10e3; 

117 Fc = 32 . 6a6 ; 

118 

119 

120 alsaif manu.sel «» 4 

121 data.f ilaname « [’ c : \resaarch\code\matlab\whites'2\wscmat"l\snr2_f 2’ wtype]; 

122 eval([’load ’ data.f ilename] ) ; 
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123 num.col = size(x,2); 

124 index.f ilename = [’c:\research\code\matlab\whites"2\cmb8_ ’ int2atr (num.col)] ; 

125 aval (['load ’ index. filename] ) ; 

126 iterations = choose(num.col,8) ; 

127 % snr = 1.6; 

128 snr * 45; 

129 Rc.DAMA * 10e3; 

130 Fc = 32 . 65e6 ; 

131 

132 

133 elseif menu.sel == 5 

134 data.f ilename = [’c:\researclAcode\matlab\whites‘2\tfscmat"l\snr2_f3’ wtype] ; 

135 aval ([’load ’ data.f ilename] ) ; 

136 num.col = size(x,2); 

137 index.f ilename = [’c:\research\code\matlab\whites-2\cmb8.’ int2str(num_col)] ; 

138 eval([’load ’ index.f ilename] ) ; 

139 iterations * choose(num_col,8) ; 

140 % snr = 1.6; 

141 snr = 45; 

142 Rc.DAMA = 10e3; 

143 Fc = 32 . 70e6; 

144 

145 

146 elseif menu.sel ■= 6 

147 data.f ilename *> [' c : \research\code\matlab\vhites"2\yscmat " l\snr2_f 5 ’ wtype]; 

148 aval ([’load ’ data.f ilename] ) ; 

149 num.col = size(x,2); 

150 index.f ilename * [ 'c : \research\code\matlab\whites"2\cmb8_ ' int2str(num.col)] ; 

151 eval([’load ’ index.f ilename] ) ; 

152 iterations = choose(num_col ,8) ; 

153 ’/.snr = 1.5; 

164 snr * 45; 

155 Rc.DAMA * 10e3; 

156 Fc = 32 . 80e6 ; 

157 

158 

159 elseif menu.sel == 7 

160 data.f ilename = [’ c: \research\code\mat lab\whites"2\wscmat"l\snr3.f 2’ wtype]; 

161 eval([’load ’ data.f ilename] ) ; 

162 num.col = size(x,2); 

163 index.f ilename = [’c:\research\code\matlab\whites"2\cmb8_’ int2str(num_ col)] ; 

164 eval([’load ’ index.f ilename] ) ; 

165 iterations = choose (num.col, 8) ; 

166 Xsnr =2.2; 

167 snr = 50; 

168 Rc.DAMA = 10e3; 

169 Fc = 32 . 65e6 ; 

170 

171 

172 else 

173 data.f ilename = [’c:\research\code\matlab\whites"2\wscmat"2\snr2.f3’ wtype]; 

174 eval([’load ’ data.f ilename] ) ; 

175 num.col = size(x,2); 

176 index.f ilename = ['c;\research\code\matlab\whites"2\cmb8_' int2str(num. .col)] ; 

177 eval([’load ’ index.f ilename] ) ; 

178 iterations = choose (num.col ,8) ; 

179 snr = 1.6; 

180 Rc.DAMA = 10e3; 

181 Fc = 32.70e6; 

182 

183 end 

184 
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185 X Some Corrections for code operation 

186 X:::::::::::::::::::::::::::::::::::::: 

187 If iterations > 1000 

188 iterations = 1000; 

189 end 

190 

191 if iterations »• 1000 

192 alpha * 1; 

193 else 

194 alpha * 2; 

195 end 

196 

197 X Load Matched Filter 

198 X;::::;:::::::::::::::: 

199 if matched.f ilter 

200 if window 

201 wtype2 - [’ha’]; 

202 else 

203 wtype2 ■ [’re’]; 

204 end 

205 

206 fname * C’pg’ int2str(Rc_DAMA/Rb_DAMA) wtype2] ; 

207 eval([’load c : \research\code\matlab\damane"l\ ’ fname]) 

208 end 


209 


210 X Process Begin 

211 X: :::::::::::::: : 

212 N « 65536; X Reset N after loading matched filter - required. 

213 X Preallocation 


214 X:::;:::::::::::: 

215 X « zeros (N, fft.avg) ; 

216 XX « zeros(N/2, fft.avg) ; 

217 X.mag.squared » zeros dengthdover.bound.index: upper .bound.index) .fft.avg) ; 

218 XX.mag. squared » zeros(N/2, fft.avg) ; 

219 X.avg » zerosdengthdower. bound.index :upper.bound_index) .fft.avg) ; 

220 XX.avg = zeros(N/2,ff t.avg) ; 

221 est.data = zeros(length(snr) ,length(lower_bound_index:upper. bound. index)) ; 

222 Xest.data ■ zerosdength(snr) ,N/2) ; 

223 

224 X Begin FFT Estimation Process 

225 X::::::::::::::::::::::::::::::: 


226 for n • 1: iterations 


227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 

244 

245 

246 


n 

temp 3 x( ; , index (n, ; )) ; 

X ■ fft(temp); 

X » Xdover.bound.index. -upper.bound. index, :) ; 
XX = Xt(l;N/2, : ) ; 

X.mag.squared » (X . *conj (X) ) . ’ ; 

X.avg = (sum(X.mag.squared) ./ fft.avg).’; 


X Select current index pattern. 
X Perform 8 FFT’s 
X Restrict to search bound 


X Sum and average . 


if matched.f ilter 

MF.mag.X * conv (MF, X.avg) ; 
end 


k » find(X_avg(l :length(X.avg)) ■» max(X_avg(l : length(X.avg) ) ) ) ; 
if insight 

dispC ’Frequency found to be:’) 

(k+20791 ) * 100e6/ 64e3 
pause(lO) 

and 

if matched.f ilter 

k * f ind(MF.mag.X(l:length(MF_mag_X)) max(MF_mag.X(l :length(MF.mag.X)))) ; 
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247 

248 if (length(k) > 1) 

249 k = median(k); 

250 end 
261 

252 k * k • Hi 

253 else 

254 k = f ind(X_avg(l:length(X.avg)) «« max(X_avg(l : length(X.avg) ) ) ) ; 

255 

266 if (length(k) > 1) 

257 k * median(k); 

258 end 

259 

260 end 

261 

262 est.dataCl ,k) = est.dataCl ,k) +1; ‘/, Increment bin 

263 

264 end 

265 

266 est_data(l , : ) = est.dataCl , : ) .» alpha; X Sponge data to look like a larger estimate 

267 

268 7. 

269 '/, Save Data in MATLAB 4.x format 

270 X 

271 x = version; 

272 if (x(l)==’5’) X running on MATLAB 5.x 

273 save -v4 cap.dat Fc Fs N N1 est.data est.tolerance number. of .estimates . . . 

274 start ing.bound snr window wtype2 

275 else X running on MATLAB 4.x 

276 save cap.dat Fc Fs N N1 est.data est.tolerance number.of .estimates . . . 

277 start ing.bound snr window wtype2 

278 end; 
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wsco_d2.m 


001 X DAMA Test Signal Analysis 

002 X 

003 X Purpose: 

004 X This code is to be used as driver code for vsco.m. This code 

005 X converts captured whites sands files and converts them into 

006 X usable length N blocks. Not entirely automated. Please read 

007 X info in vsco.m for details on project purpose. 

008 X 

009 X Input: 

010 X This code requires the DAMA test files named as sl,s2,...sn. 

011 X 

012 X Output: 

013 X Undecided at this point. 

014 X 

015 X Revision History: 

016 X 

017 X ver 1.0: baseline 

018 X ver 2.0: Changed scope of code. This program is now to 

019 X to be used to drive other code only. 

020 X ver 2.5: Restructed the format data is to be saved in. 

021 X 

022 X 

023 X Current Version - Date: 

024 X 2.5 - 8/4/98 

025 X 

026 X Author: 

027 X Brad Scaife 

028 X 

029 X Date: 

030 X 7/23/98 

031 X 

032 X Notes: 

033 X Remember to change the save filename at the end of the code. Sorry 

034 X got to do a little work to run this one. 

035 X*************** **************************** ********************************* 

036 clear all 

037 ON = 1; 

038 OFF = 0; 

039 window « OFF ; 

040 

041 X Change to working directory 

042 X****************************** 

043 cd c:\research\data\damaus-l\lk_10\snr2\f3 

044 

045 X Code Options 

046 X****************************** 

047 X capture.option » 0; X 0 for 50ksample data else 1 (for 100k) 

048 

049 X Code Parameters 

050 X****************************** 

061 num.files « 6; X 6- for snrl_f3; 9- for snr2_f3; 12- for snr3_f2 

062 half - 66536; 

053 short. cap » 60002; 

054 long.cap * 100002; 

056 

056 load si; 

057 if length(sl) =* long.cap 

058 capture.option * 1; 

059 else 

060 capture.option = 0; 
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061 and 

062 clear si; 

063 pack 

064 

065 

066 

067 

068 

069 Load raw signals 

070 %****» a***************** ******* 

071 

072 k = 0; 

073 

074 h = waitbar2(0, ‘Loading files...’); 

075 

076 while 1 

077 k ■ k+1; 

078 waitbar2(k/num_f iles) ; 

079 sk = [’s’ int2str(k)] ; 

080 filename = sk; 

081 if "exist(f ilename) , break, end 

082 eval([’load ’ filename]) 

083 end 

084 close (h) 

085 

086 '/, Preallocation 

087 %*****»• •****»***♦»* *********** 

088 i » xerosChalf , U+capture_option)*num_f iles) ; 

089 

090 

091 cd c:\research\code\matlab\vhites-2 

092 

093 if capture.option 

094 h * waitbar2(0, ’Parsing 100 kSample Files...’); 

095 index = 1; 

096 

097 if window 

098 W = hamm(half); 

099 else 

100 W * onesChalf , 1) ; 

101 end 

102 

103 for k » l:2:2*num_files 

104 waitbar2(index/num_f iles) ; 

105 varname = [‘s’ int2str C index) ] ; 

106 temp ■ eval( [varname] ) ; 

107 index * index +1; 

108 x(:,k) = temp(l .half ) . *W; 'l For others 

109 x(:,k+l) = temp(length(temp)-half+l;length(temp)) .*VI; ’/, For others 

110 end 

111 

112 close(h); 

113 else 

114 h = waitbar2(0, ’Parsing 50 kSample Files...’); 

115 

116 if window 

117 V = hammdength(sl) ) ; 

118 else 

119 W = ones(length(sl) ,1) ; 

120 end 

121 

122 for k = l:num_files 
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123 vaitbar2(k/num_f iles) ; 

124 varname * ['s’ int2str(k)] ; 

125 temp > aval ( [varname] ) ; 

126 temp = temp.*W; 

127 temp = [temp; zeros (half - length(temp) ,1)] ; 

128 x( : ,k) * temp; 

129 end 

130 close(h); 

131 end 

132 

133 cd c :\research\code\matlab\vhites'2\vscmat"2 

134 

135 mlver ■ version; 

136 if (mlver (1) =» ’5’) % Running on MATLAB 5.x 

137 save -v4 snr2_f3r x 

138 else 

139 save snr2_f3r x 

140 end 

141 

142 cd c : \research\code\matlab\vhites"2 

143 Xclear all 
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prepcap2.m 


001 X Title: 

002 X prepcap.m 

003 X 

004 X Purpose: 

005 X The purpose of this code is to prepare DAMA captured signals for 

006 % comparison against simulated results. 

007 X 

008 7. Revision: 

009 7. 1.0 — 7/26/98 

010 7. 1.5 — 8/2/98 

011 X 2.0 ~ 11/2/98 

012 X 

013 7. Revision History: 

014 X none 

015 X 1.5 — Altered code to accept neu signal matrix generated 

016 X in vsco.d.m 

017 X 2.0 — Employing Spectral Matched Filter 

018 X 

019 X Author: 

020 X Brad Scaife 

021 X 

022 X Date: 

023 X 7/26/98 

024 X 

025 X Notes : 

026 X This code requires input signals matrices that have been prepared 

027 X by vsco.d.m and vector sequence matrices prepared by goob.m. Due to 

028 X limited test set, the code makes some approximations that may or may 

029 X not be valid. Therefore, this code should not be considered in any 

030 X way a proof but rather a "ballpark" type of justification: Could it 

031 X work in the real world? 

032 X 

033 y t ************************* ****************************** *********************** 

034 clear all 
036 ON » 1; 

036 OFF - 0; 

037 

038 X Code Parameters 

039 X***************** 

040 snr = 1:7; X SNR Ranges from So/N estimate. 

041 Fc = 32.7e6; X Measured DAMA carrier - from WSC. 

042 Fs = 100e6; 7. Burr-Brown Sampling Rate. 

043 est.tolerance = 3e3; X ♦/- acceptable error. 

044 N = 65536; X FFT block length. 

045 number.of .estimates = 1000; X * of estimation attempts 

046 fft.avg = 8; 7. Number of FFT’s to perform estimation over. 

047 alpha * 2; X Correction to number of estimates. 

048 insight = 0; 

049 lower.bound. index = 21291; X Index of lover search bound. 

050 upper.bound. index = 21504; X Index of upper search bound. 

051 starting.bound = lower.bound. index; 

052 window = OFF ; 

053 matched.f ilter = ON; 

054 

056 X Preallocation 

056 X**************** 

057 X = zeros (N, fft.avg) ; 

058 

059 X.mag.squared = zeros (length(lower_bound_ index : upper.bound. index) , fft.avg) ; 

060 
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3 


061 X.avg « zeros(length(lower_bound_index:upper_bound_index) .fft.avg) ; 

062 

063 est.data = zeros (length(snr) ,length(lower_bound_index:upper_bound_index)) ; 

064 

065 if window 

066 load c:\research\code\matlab\damane'l\pglOha 

067 else 

068 load c:\research\code\matlab\damane~l\pglOre 

069 end 

070 

071 

072 X Load SNR1.F2 Data Matrix (65836x12) and Combination Matrix (495x8) 

073 %****************•**************************************•************************* 

074 load c:\research\code\matlab\uhites”2\snrl_f2 load 
076 c:\research\code\matlab\uhites~2\cmb8.12 h « 

076 waitbar2(0, ’Formulating SNR1 Estimation for f2...*); 

077 

078 for n * 1:495 

079 

080 index * cmbo3; 7. Load unique index pattern matrix. 

081 waitbar2(n/495) ; 

082 

083 temp * x( : . index(n, : )) ; X Select current index pattern. 

084 X * fft(temp); X Perform 8 FFT’s 

085 X ■ X(lower_bound_index :upper_bound_ index , : ) ; X Restrict to search bound 

086 X.mag.squared * (X.*conj (X)) . ’ ; 

087 X.avg ■ (sum(X.mag.squared) ./ fft.avg).’; X Sum and average. 

088 

089 k * find(X_avg(l:length(X_avg)) »« max(X.avg(l : length(X.avg) )) ) ; 

090 if insight 

091 disp( ’Frequency found to be:’) 

092 (k+20791 ) * 100e6/84e3 

093 pause (10) 

094 end 

096 if (length(k) > 1) 

096 k * median (k) ; 

097 disp(’ooops’) 

098 end 

099 

100 est_data(l ,k) * est.data(l,k) + 1; X Increment bin 

101 

102 end 

103 

104 est_data(l , : ) » est.datad , : ) .* alpha; X Sponge data to look like a larger estimate 

105 

106 clear x 

107 clear index; 

108 close (h); 

109 

110 

111 X Load SNR1.F3 Data Matrix (65536x12) and Combination Matrix (495x8) 

112 X************************** *********************************************** ******** 

113 load c :\research\code\matlab\whites'2\snrl_f3 

114 load c:\research\code\matlab\whites'2\cmb8.12 

115 h * waitbar2(0 , 'Formulating SNR1 Estimation for f 3 

116 

117 for n = 1:495 

118 

119 index * cmbo3; X Load unique index pattern matrix. 

120 waitbar2(n/495) ; 

121 

122 temp * x( : , index(n, : )) ; X Select current index pattern. 
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123 X = fft(temp); X Perform 8 FFT’s 

124 X = X(lower_bound_ index -.upper .bound, index, :) ; X Restrict to search bound 

125 X.mag.squared = (X.*conj(X)) . ’ ; 

126 X.avg = (sum(X.mag.squared) ./ f ft.avg) . ’ ; X Sum and average. 

127 

128 k = f ind(X_avg(l:length(X_avg)) ** max(X_avg(l:length(X_avg)))) ; 

129 if insight 

130 disp( ’Frequency found to be:’) 

131 (k+20791) *100e6/64e3 

132 pause (10) 

133 end 

134 if (length(k) > 1) 

135 k = median(k) ; 

136 disp( ’ooops ’ ) 

137 end 

138 

139 est_data(2,k) = est_data(2,k) +1; X Increment bin 

140 

141 end 

142 

143 est_data(2 , : ) - est_data(2, :) .* alpha; X Sponge data to look like a larger estimate 

144 

145 clear x 

146 clear index; 

147 close(h); 

148 

149 

150 X Load SNR2.F1 Data Matrix (65636x12) and Combination Matrix (495x8) 

151 X******************************************************************* ******** ****** 

152 load e:\research\code\matlab\uhitas"2\snr2_f 1 

153 load c:\research\code\matlab\tfhites-2\cmb8.12 

154 h « vaitbar2(0, ’Formulating SNR2 Estimation for f 1 . . . ’ ) ; 

155 

156 for n = 1:495 

157 

158 index = cmbo3; X Load unique index pattern matrix. 

159 waitbar2(n/495) ; 

160 

161 temp = x( : , index(n, : )) ; X Select current index pattern. 

162 X = fft(temp); X Perform 8 FFT’s 

163 X = X(louer_bound_index .-upper. bound.index, :) ; X Restrict to search bound 

164 X.mag.squared * (X . econj (X) ) . ’ ; 

165 X.avg « (sum(X_mag. squared) ./ fft.avg) . ’ ; X Sum and average. 

166 

167 k = find(X.avg(l:length(X_avg)) == max(X_avg(l : length(X.avg))) ) ; 

168 if insight 

169 disp( ’Frequency found to be:’) 

170 (k+20791) *100e6/64e3 

171 pause(lO) 

172 end 

173 if (length(k) > 1) 

174 k = median(k) ; 

175 disp( ’ ooops ’ ) 

176 end 

177 

178 est_data(3,k) = est_data(3,k) + 1; X Increment bin 

179 

180 end 

181 

182 est.data(3, : ) = est_data(3, :) .* alpha; X Sponge data to look like a larger estimate 

183 

184 clear x 
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185 clear index; 
188 close (h); 


187 

188 
189 


190 X Load SNR2.F2 Data Matrix (65536x12) and Combination Matrix (495x8) 


191 


••(*****«•*•*****•* 


192 load e:\research\code\matlab\vhites'2\snr2_f2 

193 load c:\research\code\matlab\whites-2\cmb8_12 

194 h • waitbar2(0, ’Formulating SNR2 Estimation for f2...’); 

195 


196 for n * 1:495 

197 

198 index » cmbo3; X Load unique index pattern matrix. 

199 waitbar2(n/49S) ; 

200 

201 temp * x index (n ; % Select current index pattern. 

202 X * fft(temp); X Perform 8 FFT’s 

203 X = X(lover_bound_index'.upper_bound_index, :) ; X Restrict to search bound 

204 X.mag.squared • (X. *conj (X) ) . ’ ; 

205 X.avg - (sum(X_mag_squared) ./ fft.avg).’; X Sum and average. 

206 

207 k « find(X_avg(l:length(X_avg)) “ max(X.avg(l : length(X.avg))) ) ; 

208 if insight 

209 disp( ’Frequency found to be:’) 

210 (k+20791)*100e6/64e3 

211 pause (10) 

212 end 

213 if ( length (k) > 1) 

214 k *> median(k); 

216 disp( ’ooops ’) 

216 end 

217 

218 est.data(4,k) > est_data(4 ,k) +1; X Increment bin 

219 

220 end 

221 

222 est_data(4, : ) = est_data(4 , : ) . * alpha; X Sponge data to look like a larger estimate 

223 


224 clear x 

225 clear index; 

226 close (h) ; 

227 

228 

229 X Load SNR2.F3 Data Matrix (65536x18) and Combination Matrix (1000x8) 

230 X* ******************************************************************************** 

231 load c:\research\code\matlab\vhites'2\snr2.f3 

232 load c : \research\code\matlab\vhites"2\mc_cmb2 

233 h « vaitbar2(0, ’Formulating SNR2 Estimation for f3...’); 

234 

236 for n * 1:1000 

236 

237 vaitbar2(n/1000) ; 

238 

239 temp * x( : , index (n, : ) ) ; 

240 X * fft(temp); X Perform 8 FFT’s 

241 X = X (lover.bound. index: upper_bound_ index ,:) ; 

242 X.mag.squared = (X . «conj (X) ) . ’ ; 

243 X.avg « (sum(X.mag.squared) ./ fft.avg).’; X Sum and average. 

244 

245 k = f ind(X.avg(l :length(X.avg)) == max(X_avg(l:length(X_avg)))) ; 

246 if (length(k) > 1) 
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247 k = median(k); 

248 disp( ’ooops ’ ) 

249 end 

250 

251 est_data(5,k) = est_data(S,k) +1; X Increment bin 

252 

253 end 

254 

255 clear x; 

256 close(h) ; 

257 

268 X Load SNR2.F5 Data Matrix (65536x12) and Combination Matrix (495x8) 

269 X******************************** ******************************************* **»*»» 

260 load c :\research\code\matlab\vhites"2\snr2_f5 

261 load c:\research\code\matlab\vhites"2\cmb8_12 

262 h = waitbar2(0, ’Formulating SNR2 Estimation for f5...’); 

263 

264 for n = 1:495 
266 

266 index = cmbo3; X Load unique index pattern matrix. 

267 vaitbar2(n/495) ; 

268 

269 temp - x( : , index(n, : ) ) ; X Select current index pattern. 

270 X * fft(temp); X Perform 8 FFT’s 

271 X = X(lover_bound_index:upper_bound_index, :) ; X Restrict to search bound 

272 X_mag_squared = (X. *conj (X)) . ’ ; 

273 X_avg = (sum(X_mag_squared) ./ fft.avg) . ’ ; X Sum and average. 

274 

275 k = find(X_avg(l:length(X_avg)) == max(X_avg(l:length(X_avg)))) ; 

276 if insight 

277 disp( ’Frequency found to be:’) 

278 (k+20791) *100e6/64e3 

279 pause (10) 

280 end 

281 if (length(k) > 1) 

282 k * median (k) ; 

283 disp( ’ooops’ ) 

284 end 

285 

286 est.data(6 ,k) * est_data(6 ,k) +1; X Increment bin 

287 

288 end 

289 

290 est_data(6, : ) ■ est_data(6, : ) .* alpha; % Sponge data to look like a larger estimate 

291 

292 clear x 

293 clear index; 

294 close(h); 

295 

296 X Load SNR3.F2 Data Matrix (65536x12) and Combination Matrix (495x8) 

297 X*********************** *************************************************** ******* 

298 load c:\research\code \matlab\vhites'2\snr3_f 2 

299 load c:\research\code\matlab\vhites"2\cmb8.12 

300 h = vaitbar2(0, ’Formulating SNR3 Estimation for f2...’); 

301 

302 for n = 1:495 

303 

304 index = cmbo3; X Load unique index pattern matrix. 

305 vaitbar2(n/495) ; 

306 

307 temp = x( : , index(n, : ) ) ; X Select current index pattern. 

308 X = fft(temp); X Perform 8 FFT’s 
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X « X(lover_bound. index :upper_bound_index ; '/■ Restrict to search bound 

A.oag.squared « (X . *conj (X ) ) . ’ ; 

X.avg = (sum(X.mag.squared) ./ fft.avg) . ’ ; X Sum and average. 


309 

310 

311 

312 

313 k » f ind(X_avg(l :length(X_avg) ) »« max(X_avg(l:length(X_avg)))) ; 

314 if insight 

315 disp( ’Frequency found to be:’) 

316 (k+20791 ) *100e6/64e3 

317 pauseClO) 

318 end 

319 if (length(k) > 1) 

320 k = median (k); 

321 disp(’ooops’) 

322 end 

323 

324 est_data(7,k) = est_data(7,k) +1; X Increment bin 

326 

326 end 

327 

328 est.data(7, : ) = est_data(7 , : ) . * alpha; % Sponge data to look like a larger estimate 

329 

330 clear x 

331 clear index; 

332 close (h); 

333 

334 

335 mlver » version; 

336 if (mlver(l) “ 6) X Running on HATLAB 5.x 

337 save -v4 dama.cap Fs Fc snr est.tolerance N number_of .estimates . . . 

338 est.data lover .bound. index upper .bound. index 

339 else X Running on MATLAB 4.x 

340 save dama.cap Fs Fc snr est.tolerance N number.of.eBtimates . . . 

341 est.data lover.bound.index upper.bound.index 

342 end 
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exl.m 


001 x********************* **************************************************** 


002 X 

003 ’/. 

004 X 

005 X 

006 % 

007 X 

008 X 

009 X 

010 */. 

011 X 

012 •/. 

013 X 

014 X 

015 */. 

016 X 

017 7. 

018 X 

019 X 

020 X 

021 X 

022 7. 

023 7. 

024 X 
026 7. 

026 7. 

027 7. 

028 7. 

029 X 

030 7. 


Experiment One: Finding Probabilities and MSE estimates of estimating 
a complex sinusoid in noise. 

Purpose: The purpose of this code is to determine and prove the 
relationship between the DAMA curves and actual MSE of 
estimation. This baseline test will provide insight into the 
simpler case of a complex sinusoid in noise which will then 
be extended to the more complex DAMA carrier case. 


Programmer: Brad Scaife 

Date: 2/14/99 

Revision Date: 4/20/99 

Current Revision: 1.10 
Revision History: 

1.0 - Baseline 

1.01 - Corrected Noise Power and theoretical curve 

1.02 - Fixed indexing problem and clearing of mse 

through each iteration. 

1.03 - Added support for estimation range compared 

to some baseline DFT resolution for comparison. 
1.10 - Increasing the computational resolution of the 

FFT so that the the sim curve can more closely 
match the Porat Curve. 


Notes: See Porat. As per indicated in Porat, the results of 

estimation are valid only when the "rule of thumb" are 
satisfied. Thus any processing of snr’s below the ROT are 
not valid with the theoretical curve. 


031 %*****•**«••****«#*************************************♦****************•* 


032 clear all 

033 clc 


034 

035 X Program Parameters : 

036 X********************* 


037 N » 512; X 

038 L = 8*N; X 

039 Nbase = N; X 

040 fs = 800e3; X 

041 Ts = l/fs; X 

042 f = 200e3 + fs/(8*N); X 

043 stime * 0; etime » (N-l)/fs; 

044 phi = 0; X 

045 A = 1; 7. 

046 Jw * 1; X 

047 D = Nbase*Ts; X 

048 

049 


050 snr = [-20: .5:30] ' ; 

051 number.estimates = 1000; 
062 


Normally set to 512. 

Computational Resolution. 

Basic DFT resolution. 

Sampling frequency. 

Sampling interval . 

Sinusoid frequency. 

Phase offset. 

Amplitude of sinusoid. 

Set as window function - see Porat. 
As per Porat . 


053 X Calculated Program Parameters 

054 X*»»***********«********»******* 


055 

056 

057 

058 


est. range = floorCf s/L)/2; 


f.lo = f - est. range; 
f.hi = f + est.range; 

059 k_lo = ceil(f_lo*N/f s) ; 

060 k_hi = f loor(f _hi*N/fs) ; 


X - Frequency estimation range as a function of 
X basic DFT resolution. 

X - Lower "accurate" estimation bound. 

X - Upper "accurate" estimation bound. 

X - Lower index bound. 

X - Upper index bound. 
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061 

062 if N *« Nbase 

063 k.hi • k.lo; 

064 and 

065 

066 res.factor = N/L; X - Resolution Factor 

067 

068 % Preallocate 

069 X************* 

070 success = zerosCl .length(snr) ) ; 

071 mse.calc * zeros ( 1 .length(snr)) ; 

072 mse * zeros (1, length (snr)) ; 

073 rmse = zeros ( 1 , length ( snr ) ) ; 

074 rmso.calc » zeros (1 , length(snr) ) ; 

075 mse.sum * 0; 

076 Xerr.sum * 0; 

077 

078 X Signal Generation, Finding of true frequency 

079 X****** **************************************** 

080 s * csin.genCf ,phi,A,fs,stime,etime) ; X Verified power * 1. 

081 signal.power = cov(s); 

082 S » fft(s) : 

083 S.mag - abs(S); 

084 Xtrue.ind * find(S_mag *» max(S_mag))-l 

085 Xf.true = true_ind/N*fs 

086 Xpause 

087 f_true « f; 

088 X true_ind = floorCf _true*N/fs) 

089 

090 X Display Parameters 

091 X******************** 

092 disp(sprintf ( ’True Frequency: XIO.Bf Hz’.f.true)) 

093 disp(sprintf ( 'Base Points: Xd', N)) 

094 disp(sprintf ( ’DFT Points: Xd’.L)) 

095 disp(sprintf ('Base DFT Resolution: Xl0.5f Hz’, fs/N)) 

096 disp(sprintf (’Calculation DFT Resolution: XIO.Bf Hz’, fs/L)) 

097 disp(sprintf (’Resolution Factor: Xd’, res.factor)) 

098 disp(sprintf( ’Estimation Range: X10.5f Hz’, est.range)) 

099 disp(sprintf (’Lower Index Bound: Xd’, k_lo)) 

100 disp(sprintf (’Upper Index Bound: Xd’, k.hi)) 

101 dispCPress a key to continue...’) pause 

102 

103 for n « 1 : length (snr) 

104 snr(n) X Current SNR to be tested. 

105 snr.mod = sqrt(cov(s)/(10"(snr(n)/10))) ; X Standard Deviation of noise 

106 X No « Ts * sqrt(2) * (snr_mod)"2; X Noise Power in W/Hz. 

107 X No « Ts * (2 * snr_mod*2)*2; 

108 No =* Ts * snr.mod'2; 

109 mse.calc(n) = (6*No*Jw)/((2*pi)"2 * A‘2 * D"3) ; 

110 X rmse.calc(n) * sqrt( 6*Jw*sqrt( cov(s)/(10‘ (snr (n)/10) ) )/(100*pi‘2))/D; 

111 X rmse.calc(n) = (l/D)*sqrt(6*Jw*PG/(100*pi'2)) ; 

112 randn( ’seed’ ,0) ; 

113 

114 for k > 1 : number. estimates 

115 v = [randn(N.l) + j *randn(N , 1)] *snr_mod/sqrt (2) ; 

116 

117 Xcov(s) X Remove me 

118 Xcov(v) X Remove me 

119 Xl0*logl0(cov(s) / cov(v)) X Remove me 

120 Xpause X Remove me 

121 y « s + v; 

122 
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123 Y = fft(y,L); 

124 Y.mag = abs(Y); 

125 XY.mag = Y_mag(l:N/2 + 1); 

126 

127 mar.ind = min (find (Y.mag ■= max(Y.mag))) - 1 

128 f .found = mar_ind/L«fs ; 

129 

130 mse.sum = mse.sum + (f.true - f_found)"2; 

131 % err.sum = err.sum + abs (f.true - f .found) ; 

132 

133 if (mar.ind >= k.lo k mar.ind <= k.hi) 

134 success(n) = success(n) +1; 

135 end 

136 

137 

138 

139 end 

140 

141 X mse(n) = (err. sum / number_estimates)"2; 

142 mse(n) = mse.sum / number.estimates ; 

143 rmse(n) = sqrt(mse(n)) ; 

144 */, err.sum = 0; 

145 mse.sum - 0; 

146 

147 success(n) - success(n)/number_estimates ; 

148 end 

149 

150 rmse.calc « sqrt(mse.calc) j 
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ex2.m 


001 %************************************************************************* 

002 X 

003 X Experiment Two: Finding Probabilities and MSE estimates of estimating 

004 X a real sinusoid in real noise. 

005 X 

006 X Purpose: The purpose of this code is to determine and prove the 

007 X relationship between the DIMA curves and actual MSE of 

008 X estimation. This baseline test will provide insight into the 

009 X simpler case of a complex sinusoid in noise which will then 

010 X be extended to the more complex DAMA carrier case. 

011 X 

012 X Programmer: Brad Scaife 

013 X Date: 3/18/99 

014 X Revision Date: 3/18/99 

015 X Current Revision: 1.0 

016 X Revision History: 

017 X 1.0 Baseline 

018 X 

019 X Notes: See Porat. As per indicated in Porat, the results of 

020 X estimation are valid only when the "rule of thumb" are 

021 X satisfied. Thus any processing of snr's below the ROT are 

022 X not valid with the theoretical curve. 

023 X********************************* ************************* *************** 

024 clear all clc 


025 


026 X Program Parameters: 

027 X********************* 


028 N « 65536: X 

029 Nbase > 65536; X 

030 fs = 800e3; X 

031 Ts » 1/fs; X 

032 f ■ 200e3 + fs/(4*N); X 

033 stime « 0; etime = (N-l)/fs; 

034 phi « 0; X 

036 A * 1; X 

036 Jw - 1; X 

037 D * Nbase*Ts; X 

038 

039 


- Normally set to 512. 

- Basic DFT resolution. 

- Sampling frequency. 

- Sampling interval. 

- Sinusoid frequency. 

- Phase offset. 

- Amplitude of sinusoid. 

- Set as window function - see Porat. 

- As per Porat . 


040 snr = [-20: .2:10] ’ ; 

041 number .estimates • 100; 

042 


043 X Calculated Program Parameters 

044 X******************************* 

045 est.range « floor (fs/Nbase)/2 ; 

046 

047 f.lo « f - est.range; 

048 f.hi « f + est.range; 

049 k.lo » ceil(f .lo*N/f s) ; 

050 k.hi « floor(f.hi»N/fs); 

051 

062 if N »» Nbase 

053 k.hi = k.lo; 

054 end 

055 

056 res.f actor = N/Nbase; 

057 

058 X Preallocate 

059 X************* 

060 success * zeros (1 , lengthCsnr)) ; 


X - Frequency estimation range as a function of 
X basic DFT resolution. 

X - Lower "accurate" estimation bound. 

X - Upper "accurate" estimation bound. 

X - Lower index bound. 

X - Upper index bound. 


X - Resolution Factor 
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061 mse.calc = zeros ( 1, length(snr) ) ; 

062 mse = zeros(l .length(snr) ) ; 

063 rase ■ zeros(l .length(snr) ) ; 

064 rase. calc = zeros(l,length(snr)) ; 

065 mse.sum = 0; 

066 Xerr.sum = 0; 

067 

068 X Signal Generation, Finding of true frequency 

069 X********************************************** 

070 s * sinu_gen(f ,phi,A,f s ,stime ,etime) ; X Verified power « 1. 

071 signal.power = cov(s); 

072 S = fft(s) j 

073 S.mag * abs(S); 

074 Xtrue.ind » find(S_mag «= max ( S.mag) )-l 

075 Xf.true » true_ind/N*f s 

076 Xpause 

077 f .true ■ f ; 

078 ‘/, true.ind = f loor (f_true»N/f s) 

079 

080 7. Display Parameters 

081 X******************** 

082 disp(sprintf (’True Frequency: 7.10 . 5f Hz’, f. true)) 

083 disp(sprintf ( ’Base Points: 7.d', Nbase)) 

084 disp(sprintf (’DFT Points: Xd’.N)) 

085 disp(sprintf ( ’Base DFT Resolution: 7.10 . 5f Hz’, fs/Nbase)) 

086 disp(sprintf ( 'Calculation DFT Resolution: Xl0.5f Hz’, fs/N)) 

087 dispCsprintf (’Resolution Factor: Xd’, res.factor)) 

088 disp(sprintf (’Estimation Range: 7,10. 5f Hz’, est.range)) 

089 disp(sprintf( 'Lower Index Bound: 7.d ’ , k.lo)) 

090 disp(sprintf( 'Upper Index Bound: Xd’ , k.hi)) 

091 disp(’Press a key to continue...’) pause 

092 

093 for n * 1 : length(snr) 

094 snr(n) X Current SNR to be tested. 

095 snr.mod * sqrt(cov(s)/(10*(snr(n)/10))) ; X Standard Deviation of noise 

096 X No * Ts * sqrt(2) * (snr.mod) "2; X Noise Power in W/Hz. 

097 X No » Ts * (2 * snr_mod"2) “2; 

098 No = Ts » snr_mod~2; 

099 mse.calc(n) = (24»No*Jw)/((2*pi)"2 * A"2 * D"3) ; 

100 7. rase.calc(n) = sqrt( 6*Jw*sqrt( cov(s)/(10‘ (snr (n)/10) ) )/ (100*pi~2) )/D ; 

101 X rase.calc(n) * (l/D)*sqrt(6*Jw*PG/(100*pi‘2)); 

102 

103 for k = 1 : number. estimates 

104 v = randn(N.l) * snr.mod; 

105 

106 X cov(s) 7, Remove me 

107 X cov(v) 7, Remove me 

108 X 10*logl0(cov(s) / cov(v)) 7. Remove me 

109 X pause X Remove me 

110 y “ s + v; 

111 

112 Y = fft(y) ; 

113 Y.mag = abs(Y); 

114 Y.mag.p = Y_mag(l:N/2 + 1); 

115 

116 max.ind = min(find( Y.mag.p == max (Y.mag.p))) - 1; 

117 f.found = max.ind/N*fs ; 

118 

119 mse.sum = mse.sum + (f.true - f_found)"2; 

120 X err.sum = err.sum + abs(f_true - f.found); 

121 

122 if (max.ind >= k.lo t max.ind <= k.hi) 
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123 suecess(n) = success(n) +1; 

124 end 

125 

126 

127 

128 end 

129 

130 % mse(n) « (err.sum / number. estimates) "2; 

131 ose(n) = mse.sum / number. estimates ; 

132 rmse(n) > sqrt(mse(n>) ; 

133 % err.sum = 0; 

134 mse.sum * 0; 

135 

136 success(n) * success Cn)/number_estimates ; 

137 end 

138 

139 rmse.calc * sqrt(mse.calc) ; 
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ex 3 .m 


001 x*******»**»**»*»*»*»*»* ****************************************** **«»*«** 

002 X 

003 X Experiment Three: Finding Probabilities and MSE estimates of estimating 

004 X a BPSK carrier in real noise. 

005 X 

006 X Purpose: The purpose of this code is to determine and prove the 

007 % relationship between the DAM! curves and actual MSE of 

008 7. estimation. This test is an extension of experiments one 

009 7. and two where a BPSK carrier is under test . 

010 7. 

011 X Programmer: Brad Scaife 

012 X Date: 3/30/99 

013 X Revision Date: 3/30/99 

014 X Current Revision: 1.0 

015 X Revision History: 

016 X 1.0 - Baseline 

017 X 

018 X Notes: See Porat. As per indicated in Porat, the results of 

019 X estimation are valid only when the "rule of thumb" are 

020 X satisfied. Thus any processing of snr’s below the ROT are 

021 X not valid with the theoretical curve. 

022 X************************************************************************* 

023 clear all clc 

024 

025 X 

026 X Program Parameters 

027 X 

028 N = 512; 

029 L « 2‘0 * N; 

030 A = 1; 

031 Jw - 1; 

032 fft.avg = 26; 

033 

034 

035 X 

036 X Communication System Parameters 

037 X 

038 Fs = 800e3; X output signal sampling freq. (samples/s) 

039 Fc = 178e3; X BPSK carrier frequency in Hz. (cycles/sec) 

040 he = floor (Fc*L/Fs) ; 

041 Rb = 10e3; X data rate (bits/s) Fs/Rb must be integer 

042 samples.per.bit = Fs/Rb; X must be integer 

043 D = N/Fs; 

044 

045 X — 

046 X Simulation Parameters 

047 X 

048 number.of. estimates « 1000;X * of frequency estimates to perform for each SNR typically 10000 

049 snr = [-12:2:14]'; 

050 msg = zeros (N,l); 

051 s * zeros(N,l) ; 

052 r » zeros(N ,1) ; 

053 

054 X — 

055 X Display Info 

056 X 

057 disp(sprintf ( 'Sampling Frequency Fs : X5. 15f ’ ,Fs) ) ; 

058 disp(sprintf ( 'Carrier Frequency Fc: X5 . 15f ' ,Fc) ) ; 

059 disp(sprintf CFFT Resolution (Data Supported): 7.5 . 15f ’ ,Fs/N)) ; 

060 disp(sprintf ( 'FFT Computational Resolution: 7,5. 15f ’ ,Fs/L)) ; 
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061 

062 

063 pause 

064 

065 X— — 

066 X Begin Iterative SNR Loop 

067 X 

068 

069 for k = 1 : length(snr) 

070 

071 rand( ' seed’ , 1000) ; 

072 randn( ’seed’ ,0) ; 

073 mse.sum » 0; 

074 

075 noise.pouer * (A"2/2)/(10"(snr(k)/10)) ; 

076 No = (noise.pouer) / Fs; 

077 help_f actor = 1; 

078 mse_calc(k) » (24*No«Jv)/(help_factor*(2*pi)‘2 * (A"2/(4*Rb)) * D"3) ; 

079 

080 for 1 = 1 :number_of_estimates 

081 

082 R.mag.sum = zeros(L/2, 1) ; 

083 

084 for n = l.’fft.avg 

085 X 

086 X Generate BPSK Signal 

087 X 

088 msg ■ A*[eos(2*pi*( [0:N-1] )*Fc/Fs)] . ’ ; 

089 data = filter(ones(samples_per_bit,l) ,1 .upsamp ... 

090 ( (-1) . "(round(rand(ceil{N/samples_per_bit) , 1))) , . . . 

091 samples.per.bit)) ; 

092 

093 data « data(l :length(msg) ) ; 

094 s » msg .* data; 

095 

096 

097 noise * randn(length(msg) ,1) ,* sqrt(cov(s)/(10"(snr(k)/10)) ) ; 

098 r = s + noise; 

099 r = r ./ sqrt(cov(r)) ; 

100 
101 
102 

103 R « fft(r,L); 

104 R_mag_sum = R.mag.sum + R(l:L/2) .*conj(R(l:L/2)) ; 

105 end 

106 

107 R_mag * R_mag_sum ./ fft.avg; 

108 

109 kmax » min(f ind(R.mag »■ max(R.mag) ) ) ; 

110 f_ast ■ (kmax-l)/L*Fs ; 

111 

112 mse_sum = mse_sum ♦ (Fc - f_est)"2; 

113 XCFc - f _est p kmai] 

114 

115 end 

116 

117 mse(k) =» mse_sum/number_of .estimates ; 

118 rmse(k) « sqrt(mse(k) ) ; 

119 rmse_calc(k) * sqrt (mse.calc(k) ) ; 

120 end 

121 

122 elf 
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123 plot(»nr,rmse_calc, ‘ — ’) 

124 hold on 

126 plot(snr,rmse) 

126 hold off 

127 grid 
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B Motorola DSP 56303EVM Code 
rev30.asm 


001 ; REV 3.0 

002 

003 ; Just lor convsniencs - delete later! !! Turns on/off D/A codec 

004 ; bclr #19,x:M_CRB0 ;Disable Rx on A-codec 

005 ; bclr #18,x:M_CRB0 ;Used to disable Tx interrupt 

006 

007 

008 opt now 

009 

010 ; nolist 

011 include 'ioequ.asm' 

012 include 'intequ.asm' 

013 include ’ada.equ.asm’ 

014 include 'vectors .asm’ 

015 include ’7819equ.asm’ 

016 ; list 

017 

018 ; ****************************************************************************** 

019 ; Initial Layouts: This section of code sets up the D/A memory resources, 

020 ; the program memory resources and defines the FFT macro. 

021 ; ****************************************************************************** 


022 

include 

’CS4216.asm' 

;D/A Memory Resources 


023 

include 

’fftr2cn.asm’ 

;FFT Macro 


024 

include 

’convm.asm’ 

; Convolution Macro 

TT - 

025 

026 

include 

'mlayout.dat' 

-.Memory Layout 



027 ; ****** ************ ***************************************** ****** ************* 

028 ; Fast Interrupt - IRQB 

029 ;****************************************************************************** 

030 org pli:I_IRQB 

031 movep y:BB7819_DR,x: (r0)+ 

032 org pli : I.IRQB+l 

033 bset #0,x: FLAGS 

034 

035 org p:$100 

036 START 

037 main 

038 

039 ; ****************************************************************************** 

040 ; Set Operating Frequency 

041 ;****************** ************ ************************ ************* *********** 

042 movep #CLK_RATE,x:M_PCTL ;Set PLL and Chip Operating Frequency 

043 

044 ;**********•************ ************** ************ *************** ************** 

045 ; Set Operating Parameters of DSP56303 

046 ; ******* *•••••*•*•****•*****•***••*****•*»*•*•******•**•*•******• *••*******•*•• 

047 move *0P_M0DE,omr ;Set Operating Mode of 303 

048 

049 ; ****************************************************************************** 

050 ; Setup Stack 

051 ;****************************************************************************** 


052 

movec 

*0,sp 

; clear hardware stack pointer 

063 

move 

•STACK, r6 

; initialise stack pointer 

054 

move 

*-l,m6 

;linear addressing 

055 




056 




057 





058 ;****************************************************************************** 
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059 ; Set AAR Wait States for External Memory(32k) and A/D Codec 

060 ;****************«»****** ********* ************************* ••*******••*• ******* 

061 include ’vs.set.asm’ ;Set Wait States 

062 

063 ; ****************************************************************************** 

064 ; Set Up IRQB Interrupt Parameters. IRQB is the interrupt designated to 

065 ; the A/D Codec. 

066 ;*****#*****#*******#**#**#**************************#************************* 


067 

ori 

#$03, mr 

; Mask all interrupts until needed 

068 

include 

’ core. ipl. asm’ ; 

i Set IRQB Interrupt Parameters 

069 




070 




071 

AAR2 equ 

$fff c21 ; 

; Compare Upper 12 bits to fffxxx 

072 

movep 

#AAR2,x:M_AAR2 ; 

; Setup AAR2 

073 




074 




075 

jsr 

INIT ; 

; Register Initialization Routine 

076 

andi 

#$fc,mr ; 

; Re-enable all interrupts 

077 




078 

main.loop 



079 




080 

jclr 

#0 , x : FL AGS , * ; 

;Wait for Sample In 

081 

bclr 

#0,x: FLAGS 


082 

jsr 

process.sample 


083 

jmp 

main.loop 


084 




085 

process.sample 



086 

clr 

a 


087 

move 

x : SAMPLE.COUNTER 

;,a0 

088 

dec 

a 


089 

move 

aO , x : SAMPLE.COUNTER 

090 

tst 

a 


091 

jne 

GET.NEXT.SAMPLE 


092 




093 

ori 

#$03, mr ; 

Disable Interrupts 

094 




095 

jsr 

WIN.N.SCALE 


096 

jsr 

COMPUTE.FFT 


097 

jsr 

AVG.FFT 


098 




099 

andi 

#$fe,mr 

; Enables only 800k A/D 

100 

movep 

#IPRC ,x :M_IPRC 

; Re-enable A/D Codec 

101 




102 

GET.NEXT.SAMPLE 



103 

rts 



104 




108 

; Subroutines: 



106 

;**••***•**•*** 



107 

include 

‘ comp.f f t . asm’ 


108 

include 

’wsc .asm' 


109 

include 

’avg.fft .asm’ 


110 

include 

’init.asm’ 


111 

include 

’get.bin. asm’ 


112 

include 

’ sinvgid . asm’ 


113 

include 

’ada. init.asm’ 


114 

echo 



115 

end 
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mlayout.dat 


001 ;******************** ****** *************** *************************** ********** 

002 {mlayout.DAT: This data file is used vith revl.ASM to lay things out in memory 

003 ; 

004 ; Notes: For use vith rev 2.1 code. 

005 ;****************************************************************************** 

006 ; References: 

007 ; DSP66300 Family Manual (300FM) 

008 ; DSPB6303 User’s Hanual (303UM) 

009 ; ****************************************************************************** 


010 ; Equates: 

011 ; ********** 

012 

013 ; ****************************************************************************** 


014 ; 56303 Processor Operating Parameters Control 

015 ; **•**•***•••*****•***••*•****•••*******•*****•••*****•****•**••*•***•*•••*••*• 

016 CLK.RATE equ $040004 ;Chip Operating Clock - See Section 9.3 300FM 

017 OP.MODE equ $389 ;Chip Operating Mode 

018 ; Please use either 389 or 3C9 for proper 

019 '.operation. Please see DAMA Programming Notes 

020 ;and 303UM:3-13 for details. 

021 FS equ 32000 ;Please set the same as D/A sample rate. 

022 ; •****•*•*••••*••***********•••••**•**••***••**•**•**•**••***•••******•***•••** 


023 ; DAMA Project Memory Settings. PLEASE DO NOT CHANGE!!! CODE WILL LIKELY NOT 

024 ; FUNCTION. THE HEMORY HAS BEEN SPECIFICALLY SETUP UTILIZING ALL ON-CHIP 

025 ; MEMORY. 

026 ; ****************************************************************************** 


027 

POINTS 

®qu 

512 

-.Number of Points (samples) 

028 

TABLE.SIZE 

equ 

612 

;Sine Wave Lookup Table Size (Will adjust output) 

029 

ITERS 

equ 

8 

; FFT Iterations 

030 

OFFSET 

equ 

128 

; Correction from Spectral Smearing due to Convolut 

031 

;OUTPUT_SEC 

equ 

2 

; Please enter duration of output in seconds. 

032 

{OUT.LENGTH 

equ 

CCVI (FS*OUTPUT_SEC) 

033 





034 

; Long Memory: 




035 

;************** 




036 

org l:$OOOa 




037 





038 

SAMPLE.DATA 

dsm 

POINTS ; Signal buffer (0200 - 03ff) 

039 

FFT.DATA 

dsm 

POINTS ;FFT Output buffer (0400 - 05ff) 

040 

;FFT .RESULT 

dsm 

POINTS ; Result FFT Data (0600 - 07ff) 

041 

COEFF 

dsm 

POINTS jSine-Cosine "Twiddle" Factor Lookup (0800 - 09ff) 

042 





043 

; X Memory : 




044 

; *********** 




045 

org i:$000a 

; see 

ADA.INIT.ASM for why we start at x:$000a 

046 





047 

SA.DATA.PTR 

ds 

1 

jSAMPLE.DATA Pointer Storage 

048 

FT.DATA.PTR 

ds 

1 

jFFT.DATA Pointer Storage 

049 

IFFT.PTR 

ds 

1 

; Imaginary FFT Data Pointer Storage 

050 

IFFT.MOD 

ds 

1 

; Imaginary FFT Data mod Storage 

051 

jFT.RES.PTR 

ds 

1 

jFFT.RESULT Pointer Storage 

052 

MAG.PTR 

ds 

1 

{Magnitude Squared Data Pointer Storage 

053 

COEFF.PTR 

ds 

1 

jCoeff Pointer Storage 

064 

WAV.PTR 

ds 

1 

jSine Wave Table Pointer Storage 

055 

WIN.PTR 

ds 

1 

jWindov Pointer Storage 

056 

WIN. MOD 

ds 

1 

{Window Modulo Storage 

057 

SMF.PTR 

ds 

1 

;SMF Pointer Storage 

058 

CNVO.PTR 

ds 

1 

.'Convolution Result Buffer Pointer Storage 

059 

FFT.COUNTER 

ds 

1 

;FFT Counter 

060 

SAMPLE. COUNTER 

ds 

1 

{Sample Counter 
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061 

MAX.VAL 

ds 

i 

Maximum value storage 

062 

MAX.LOCATION 

ds 

i 

Holds address of max location 

063 

INT.DELTA 

ds 

i 

Delta for Carrier Reconstruction 

064 

FRAC .DELTA 

ds 

i 



065 

RO.STORE 

ds 

i 

rO 

storage 

066 

Rl. STORE 

ds 

i 

rl 

storage 

067 

R2.ST0RE 

ds 

i 

r2 

storage 

068 

R3.ST0RE 

ds 

i 

r3 

storage 

069 

R4.ST0RE 

ds 

i 

r4 

storage 

070 

R5.ST0RE 

ds 

i 

r5 

storage 

071 

R7.ST0RE 

ds 

i 

r7 

storage 

072 

MO. STORE 

ds 

i 

mO 

storage 

073 

M1.ST0RE 

ds 

i 

ml 

storage 

074 

M2.ST0RE 

ds 

i 

m2 

storage 

076 

M3.ST0RE 

ds 

i 

m3 

storage 

076 

M4.ST0RE 

ds 

i 

m4 

storage 

077 

M5.ST0RE 

ds 

i 

m5 

storage 

078 

M7_ STORE 

ds 

i 

m7 

storage 

079 

N5.ST0RE 

ds 

i 

n5 

storage 

080 

N7.ST0RE 

ds 

i 

n7 

storage 

081 

FLAGS 

ds 

i 

User Defined Flag Register 

082 

OUT.COUNTER 

ds 

i 

Output Sample Counter 

083 

CNV.MEM 

ds 

i 

For use in convolutional code. 

084 

STACK 

equ 

* 

Beginning of Stack 


085 

086 org x : $800 

087 ;**************♦*»**************** 

088 ; Magnitude Squared Data 

089 ; ********************************* 

090 MAG.Sq.DATA dsm POINTS/2 

091 

092 org x:$A00 

093 ;**********»*****»*«****«********» 

094 ; Generate Sine Wave Lookup Table 

095 ; ********************************* 

096 TAB dsm TABLE.SIZE 

097 include 'sintab.asm' 

098 s intab T ABLE.SIZE , T AB 

099 

100 ; Y Memory: 

101 ; *********** 

102 org y:$0 

103 ; ************************************* 

104 ; Generate Hamming Window v/ Prescale 

105 ; ************************************* 

106 HAMM dsm POINTS ; Hamming Window table. 

107 include 'hamming. asm' 

108 hamming POINTS ,HAMM 

109 

110 

111 ; Build Twiddle factor lookup tables for FFT Routine 

112 ;***************** *************************** ******** 

113 include ’sincos.asm’ [Twiddle factor macro - builds lookup tables 

114 sincos POINTS, COEFF [Build lookup tables. 

115 

116 

117 ; Spectral Matched Filter 

118 ;***»********************* 

119 SMF dsm POINTS/2 

120 org y:SMF 

121 include 'smf20.dat’ 

122 
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123 

124 ; Convolution Output 

125 ;**•»**»•*****«*»**** 

126 CNV.OUT dam 2*P0INTS-1 



wsc.asm 


01 

02 

03 

04 

05 

06 

07 

08 

09 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 


win_n_scale Subroutine 

Purpose: The purpose of this subroutine is to scale the input data to avoid 
overflow problems. This will have the effect of lowering the overall values 
of the spectrum but will not alter the shape of the spectrum. 

In: x:(r0) - Sample Buffer 

y:(r4) - Hamming buffer 

OUT: x:(r0) - Sample Buffer w/ Window, scale, and iteration 

adjustments 

Alters: a,b,r0,r3,r4,x0,xl,y0,yl 

Written By: Brad Scaife 

Date: 2/20/98 

Platform: Motorola DSP56303 

Calls: None 

This code is verified for use with version three code. See rev30.asm. 
****************************************************************************** 
WIN.N. SCALE 


move x:(r0)+,x0 y:(r4)+,y0 ;Preload values, 

do #P0INTS/2-l,ND_SCALE ; 

mpyr xO,yO,a x:(r0)+,xl y:(r4)+,yl ;x’ (x)*(w(n)*scale/# iterations 

mpyr xl,yl,b x:(rO)+,xO y:(r4)+,y0 jSecond Iteration 

move a,x:(r3)+ ;Store a into sample buffer 

move b,x: (r3)+ , -Store b into sample buffer 

ND.SCALE 

mpyr xO,yO,a x:(rO)+,xl y:(r4)+,yl ;Loop clean up: Two mults and 
mpyr xl,yl,b ; corresponding writes to memory 

move a,x: (r3)+ ;Counters back to top of 

move b,x:(r3)+ ibuffer 
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compJft.asm 


01 ;*******•********************************* **************************+********** 

02 ; COMPUTE FFT Subroutine 

03 i 

04 ; Purpose: The purpose of this subroutine is to compute the FFT of the input 

05 ; signal and store it in memory. 

06 ; 

07 ; In: r0,rl,r2,r3,r4,r6,m0,ml ,a2,m3,m4,m5 

08 ; Out: x: (rl) ) rO,rl I r2,r3,r4 > rS,mO,ml ,m2 l m3,m4,m5 

09 ; Alters : Everything 

10 ; ****************************************************************************** 
11 COMPUTE.FFT 


12 

move 

rO , x : RO.STORE 

13 

move 

rl,x:Rl_ STORE 

14 

move 

r2 ,x :R2. STORE 

15 

move 

r3,x:R3_ST0RE 

16 

move 

r4 , x : R4_ST0RE 

17 

move 

r5,x:R5.ST0RE 

18 

move 

r 7 , x : R7_ST0RE 

19 

move 

mO , x : MO. STORE 

20 

move 

ml ,x :M1_ST0RE 

21 

move 

m2,x:M2_ST0RE 

22 

move 

m3, x: M3 .STORE 

23 

move 

m4 , x : M4.ST0RE 

24 

move 

mS , x : M5.ST0RE 

25 

move 

m7 , x : M7.ST0RE 

26 

move 

n5,x:N5. STORE 

27 

move 

n7,x:N7. STORE 


28 

29 ff tr2cn POINTS , SAMPLE.DATA ,FFT_DATA , COEFF 

30 


31 

move 

X :N7_ST0RE,n7 

32 

move 

X : N5.ST0RE , n5 

33 

move 

x:M7_ST0RE,m7 

34 

move 

x:H5_ST0RE,m5 

35 

move 

x :M4. STORE ,m4 

36 

move 

x:M3_ST0RE,m3 

37 

move 

x: M2. STORE, m2 

38 

move 

x: Ml. STORE, ml 

39 

move 

X :MO.STORE,mO 

40 

move 

x:R7_ST0RE,r7 

41 

move 

x:R5.ST0RE,r5 

42 

move 

x:R4_ST0RE,r4 

43 

move 

x:R3_ST0RE,r3 

44 

move 

x :R2_ STORE, r2 

45 

move 

x :Rl_STORE,rl 

46 

move 

x: RO.STORE, rO 

47 

rts 
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avgjft.asm 


01 ; for use vith rev 2.1 code only AVG.FFT 

02 


03 


move 

x : (rl)+ ,x0 y : (r7)+,y0 

04 


do 

#(P0INTS/2) .END.TLOOP 

OS 


mpy 

xO.xO.a x : (r2) ,yl 


06 


macr 

yO.yO.a 


07 


add 

yl.a 


08 


move 

x: (rl)+,xO y: (r7)+,y0 

09 


move 

a,x:(r2)+ 


10 





11 

END.TLOOP 



12 


clr 

b 


13 


move 

x : FFT_C0UNTER , bO 


14 


dec 

b 


15 


move 

bO , x : FFT.CQUNTER 


16 


tst 

b 


17 


jseq 

GET_MAX.BIN 


18 





19 

• **************************************** 

20 

; Prepare to 

perform next FFT iteration 

21 

• **************************************** 

22 


move 

•POINTS, xl ;Reload sample counter for next sample 

23 


move 

xl , x : SAMPLE_COUNTER 

; buffering. 

24 





25 

* 

move 

#FFT_RESULT,r2 

; re-Setup FFT Result ptr 

26 


move 

#MAG_SQ_DATA,r2 

; re-Setup Mag Squared data ptr 

27 


move 

*FFT_DATA,rl 

; Setup FFT Data ptr 

28 


move 

rl,r7 

;Imag. Pointer to FFT Buffer 

29 





30 


move 

#*0,x0 


31 


do 

fPOINTS.CLR.DAT 

'.Clear FFT Data buffer 

32 


move 

x0.x:(rl) 

; Real 

33 


move 

xO,y : (rl)+ 

; Imaginary 

34 

CLR_DAT 




35 





36 


do 

•POINTS , CLR_SHP 


37 


move 

xO,x: (rO) 


38 


move 

xO,y : (r0)+ 


39 

CLR.SMP 




40 


rts 
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get.bin. asm 

01 •, ****************************************************************************** 

02 ; GET_MAX_BIN Subrout ine 

03 i 

04 ; Purpose: The purpose of this subroutine is to determine the frequency bin 

05 ; that has the largest component and then to determine the delta for the 

06 ; sine wave generation routine. 

07 ; 

08 ; In: x:(r2) 

09 : Out : b 

10 ; Alters: b,xl,yl > r2 

11 i 

12 ; Notes: For use with rev 2.1 code only!! 

13 ;*************•*+*•*****»**************** ******* ******************************* 

14 GET-MAX_BIN 

15 

16 ; **************************************** 

17 ; Clean Up From GET-BIN Subroutine 

lg ; **************************************** 


19 

move 

iPQINTS'Xl ; Reload sample counter for next 

20 

move 

xl.x: SAMPLE-COUNTER jbuffering. 

21 

move 

#FFT_DATA,rl ; re-Setup FFT Data ptr 

22 

move 

rl,r7 ;Imag. Pointer to FFT Buffer 

23 

; move 

*FFT_RESULT , r2 ;re-Setup FFT Result ptr 

24 

move 

tMAG.SQ_DATA.r2 ;re-Setup Mag Sq Data ptr 

25 

J ************4 


26 

; Perform SMF 

Convolution 

27 

;****************************** 

28 

move 

rO,x:RO. STORE 

29 

move 

rl ,x:Rl. STORE 

30 

move 

r4 , X : R4_ST0RE 

31 

move 

mO,x: MO-STORE 

32 

move 

ml, x: M1_ST0RE 

33 

34 

move 

m4 , x : M4-ST0RE 

35 

36 

convm 

POINTS/2-1, MAG-SQ-DATA , SMF , CNV.OUT , CNV.HEM 

37 

move 

x : RO-STORE , rO 

38 

move 

x:Rl_STORE,rl 

39 

move 

x : R4-ST0RE , r 4 

40 

move 

x : MO .STORE ,m0 

41 

move 

x :Ml.STORE,ml 

42 

43 

move 

x:M4.ST0RE,m4 

44 

move 

•CNV-OUT , r7 

45 

move 

*P0INTS/2-2,m7 

46 

clr 

b 

47 

do 

•POINTS-1, ND.MAX 

48 

move 

x:(r7),xl ;Bin Comparison 

49 

cmp 

xl,b ib-xl 

50 

jit 

NEW-MAX ;b will always hold max 

51 

52 

jmp 

DUMMY 

53 

NEW-MAX 


54 

move 

xl,b ;Store New Max Location 

55 

move 

r7 , x : MAX-LOCATION 

56 

DUMMY 


57 

move 

(r7) + 

58 

nop 


59 

ND.MAX 


60 

move 

x: MAX-LOCATION, b ;Subtract mu location from 
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123 ; 




124 

NOP 


Reserved 

125 

NOP 



126 ; 




127 

j®P 

* 


128 

NOP 


SCI Receive Data 

129 ; 




130 

jmp 

* 


131 

NOP 


j- SCI Receive Data w/ Except 

132 ; 




133 

jmp 

* 


134 

NOP 


SCI Transmit Data 

135 ; 




136 

jmp 

♦ 


137 

NOP 


SCI Idle Line 

138 ; 




139 

jmp 

* 


140 

NOP 


SCI Timer 

141 ; 




142 

NOP 


Reserved 

143 

NOP 



144 | 




145 

NOP 


Reserved 

146 

NOP 



147 ; 




148 

NOP 


Reserved 

149 

NOP 



150 ; 




151 ; 




152 

j“P 

* 


153 

NOP 


; Host receive data full 

154 ; 




155 ; 




166 

jmp 

* 


157 

NOP 


Host transmit data empty 

158 ; 




159 

jmp 

* 


160 

NOP 


; Available for Host Command 

161 

jmp 

* 


162 

NOP 


; Available for Host Command 

163 

jmp 

* 


164 

NOP 


; Available for Host Command 

165 

jmp 

* 


166 

NOP 


i Available for Host Command 

167 

jmp 

* 


168 

NOP 


; Available for Host Command 

169 

jmp 

* 


170 

NOP 


; Available for Host Command 

171 

j“P 

* 


172 

NOP 


; Available for Host Command 

173 

jmp 

* 


174 

NOP 


; Available for Host Command 

175 

jmp 

* 


176 

NOP 


; Available for Host Command 

177 

jmp 

* 


178 

NOP 


; Available for Host Command 

179 

jmp 

* 


180 

NOP 


; Available for Host Command 

181 

j“P 

* 


182 

NOP 


; Available for Host Command 

183 

jmp 

* 


184 

NOP 


; Available for Host Command 
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61 ; 

move 

tFFT_RESULT-l ,yl 

62 

move 

#MAG_SQ_DATA-1 ,yl 

63 

sub 

yl.b 

64 

move 

tOFFSET , y 0 

66 

sub 

yO.b 

66 

move 

bl ,nS 

67 

j*r 

SIHWGID 

68 



69 

move 

#>ITERS,xl 

70 

move 

x 1 , x : FFT.COUNTER 

71 

rts 



location to get the actual index 
location to get the actual index 
Equals index imax 


;end program 
;Init FFT Counter 

iReset FFt counter for next iteration 
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sinwgid.asm 


01 

**** *************** ************************** ******* ********** **************** 

02 

SINWGID Subroutine 


03 




04 

Purpose: The 

purpose of this subroutine is to generate a tone at a frequency 

05 

based upon the delta value passed 

in from GET. BIN. ASM. 

06 




07 

In: b 



08 

Out: n/a 



09 

Alters: a,xl,r5,n5,y0, 


10 




11 

Notes : 



12 

For use with rev 2.1 code. 


13 

********* **** *************** ************************************************** 

14 SINWGID 



15 




16 

movep 

#$0,x:M_IPRC 

;Disable A/D 

17 

andi 

#$fc,mr ; Enable all interrupts 

18 




19 

****************************************************************************** 

20 

Initialization: D/A Codec and Setup Control Words. Only initialize the 

21 

first time. 



22 

****************************************************************************** 

23 

jset 

#l,x:FLAGS,send_loop 

;Skip after initial pass 

24 

jsr 

ada.init 

; initialize codec 

25 




26 

move 

#T0NE_0UTPUT,y0 

;set up control words 

27 

move 

y0,x:TX_BUFF_BASE+2 


28 

move 

•TONE. INPUT, yO 


29 

move 

yO , x : TX.BUFF.BASE+3 


30 




31 

bset 

#l,x: FLAGS 

;Set after initialization. 

32 




33 




34 send.loop 



35 

jset 

#2 »x:M_SSISR0 ;vait for frame sync to pass 

36 

jclr 

#2,x:M_SSISR0,* ;vait for frame sync 

37 




38 

move 

x: (r5)+n5,y0 


39 




40 

clr 

a 

;Test for end of duration of 

41 

move 

x : 0UT.C0UNTER , aO 

samples out phase. 

42 

dec 

a 


43 

move 

aO , x : 0UT.C0UNTER 


44 

tst 

a 


45 

j®q 

restart 


46 




47 

move 

yO , x : TX.BUFF.BASE 

transmit left 

48 

move 

yO,x:TX.BUFF.BASE+l 

transmit right 

49 

j®P 

send. loop 


50 




51 restart 



52 

move 

#0UT.LENGTH,xl 

ReSet output duration counter. 

53 

move 

x 1 , x : 0UT.C0UNTER 


54 




55 

rts 
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wsjset.asm 


01 

02 

03 

04 
06 
06 

07 

08 

09 

10 
11 
12 
13 


•* a******************************************* ********* **»«***»**«»*»*****♦**» 
Wait State Parameter Settings Routine 

Purpose: The purpose of this routine is to set the Bus Control Register (BCR) 
to the proper number of vait states required by each AAR device. In the 
DAMA Project, the AAR devices are: 

32k SRAM 

Burr-Brown Codec (Operating # 800kHz for DAMA Project) 

****************************************************************************** 
Wait State Settings : 

♦a****************** 


14 

DEFAULT _WS 

equ 

$0f 

■.default are vait states 

15 

SRAM.WS 

equ 

♦Of 

;32KW SRAM 

16 

FLASH.WS 

equ 

♦00 

; FLASH 

17 

PERIPH.WS 

equ 

♦Of 

;A/D Peripheral board 

18 





19 





to 

o 

AREA0 

equ 

SRAM.WS 


21 

AREA1 

equ 

FLASH. WS 


22 

AREA2 

equ 

PERIPH.WS 


23 

AREA3 

equ 

SRAM.WS 


24 





26 

BBS 

equ 

♦0 

;Bus State 

26 

BLH 

equ 

♦0 

: Bus Lock Hold 

27 

BRH 

equ 

♦0 

.Bus Request Hold 


28 

29 

30 BCR 

31 

32 

33 


equ (BBS«21)+(BLH«22)+(BRH«23)+(DEFAULT_WS«16kM_BDFW)\ 

+ (AREA3<<13kH_BA3W) + (AREA2<<10AM_BA2W) + (AREAl<<5fcM_BAlW) + (AREAOItM_BAOW) 

movep #BCR , x : M_BCR ; Initialize Bus Control Register 
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convm.asm 


convm macro length, xcoefs .hcoefs .result ,cnv_mem 

convm ident 1,0 

;****************************************«**** 4^***jM****** ************* 

t 

; Macro Name: CONVM.ASM 


Purpose: The purpose of this macro is to provide the convolution 

of two sequences stored in memory. The algorithm does 
a nested structure to minimize the memory required. 


Programmer: 
Initial Date: 
Current Rev: 

Curr. Rev. Date: 
Revision History: 


Brad Scaife 
2/20/99 
1.0 

2/20/99 


1.0 - Baseline 


Legal Statement : 

This DSPS6xxx macro may be freely used with out the permission 
of the author. The author provides the code vith the intent that 
it is not to be used where such use may endanger life and property. 
Use of this macro code releases the author from ANY litagation both 
past, present, and future from ANY and ALL such liability claims. 
Use of this code is expressly permitted at your own risk. 


Resources Used: 

Registers Used: 

a,b,r0.rl,r4,n0,n4,m0,ml,m4,x0,xl,yl 


Notes: 

Please note that this revision of the code requires the two 
input sequences to be of equal length. 

44444444444444444444444444444444444444444444444444444444444444444444444 


K equ length 

K.ALL equ 2*K-1 


move 

move 

move 

move 

move 

move 

; Begin Calculation 

move 

clr 

clr 

move 

do 

move 

do 


txcoef s ,r0 
#K-l.m0 
thcoef s ,r4 
#K-l,m4 
•result ,rl 
•K. ALL, ml 


•0,x0 

b x0,x:CNV_MEM 


x:(r0),x0 y:(r4),yl 
•K.ALL/2+1, FIRST 
b0,x : CNV.MEM 
x:CNV_MEM,END_F 
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END.F 


mac xO,yl,a x:(rO)-,xO y:(r4)+,yl 


inc b 

move <hcoefs,r4 

macr xO.yl.a 

move bO , nO 

move a,x:(rl)+ 

move txcoefs,rO 

nop 
nop 

clr a 

move CrO)+nO 

move x:(rO)-,xO y:(r4)+,yl 

FIRST 

dec b 

dec b 

move #>l,xl 

move x:(rO)-,xO y:(r4)+,yl 

do #K_ ALL/2, LAST 

move bO , x : CNV.MEM 

do X : CNV.MEM , END _L 

mac xO,yl,a x:(rO)-,xO y:(r4)+,yl 

END.L 

move #hcoefs+l,r4 

move xl,n4 

move b0,x:(r6)+ 

move xl.bO 

inc b 

macr xO,yl,a 

move bO,xl 

move #xcoefe+K-l,rO 

move x:-(r6),b0 

move a,x:(rl)+ 

dec b 

clr a 

move (r4)+n4 

move x:(rO)-,xO y:(r4)+,yl 


LAST 

endm 
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cs4215.asm 


01 ; Buffer for talking to the CS4215 

02 


03 

org 

x : 0 









04 

RX. BUFF. BASE 

equ 

* 








05 

RX.data. 1.2 

ds 

i 

,'data 

time 

slot 

1/2 

lor 

RX 

ISR 

06 

RX.data.3.4 

ds 

i 

;data 

time 

slot 

3/4 

for 

RX 

I SR 

07 

RX_data_5.6 

ds 

i 

; data 

time 

slot 

5/6 

for 

RX 

ISR 

08 

RX.data_7_8 

ds 

i 

,-data 

time 

slot 

7/8 

for 

RX 

ISR 

09 











10 

TX.BUFF.BASE 

equ 

* 








11 

TX.data.1.2 

ds 

i 

;data 

time 

slot 

1/2 

for 

TX 

ISR 

12 

TX_data_3_4 

ds 

i 

;data 

time 

slot 

3/4 

for 

TX 

ISR 

13 

TX_data_5_6 

ds 

i 

;data 

time 

slot 

5/6 

for 

TX 

ISR 

14 

TX_data.7_8 

ds 

i 

;data 

time 

slot 

7/8 

for 

TX 

ISR 


15 

16 RX_PTR 

17 TX.PTR 

18 

19 T0NE.0UTPUT 

20 TONE. INPUT 

21 CTRL.WD.12 

22 CTRL.VfD.34 

23 CTRL.WD.56 

24 CTRL.VfD.78 


ds 

ds 

EQU 

EQU 

equ 

equ 

equ 

equ 


Pointer for rx buffer 
Pointer for tx buffer 


HEADPHONE.EN+LINEOUT_EN+(4*LEFT_ATTN)+(4*RIGHT_ATTN) 

MIC_IN_SELECT+(15*MONITOR_ATTN) 

NO.PREAMP+HI.PASS.FILT+SAMP.RATE.32+STEREO+DATA.16 

IMMED.3STATE+XTAL1.SELECT+BITS.64+C0DEC.MASTER 

$000000 

$000000 


;CLB=0 
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init.asm 


01 ;**»**» ****************************************** ****************************** 

02 ; INIT Subroutine 

03 ; 

04 ; Purpose: The purpose of this subroutine is to initialize pointers to memory 

05 ; and clear out buffers . 

06 ; 

07 ; In: none 

08 : OUT: none 

09 ; Alters: r0,m0,rl,mi,r2 t m2,r3,m3,r4,m4 

10 ; 

11 ; Notes: For use with rev 3.0 code. 

12 ;******************»**************************» ******************************** 
13 INIT 


14 

move 

#S AMPLE_DATA.ro 

Setup Sample buffer ptr 

15 

move 

r0,r3 

Alternate Sample buffer ptr 

16 

move 

#P0INTS-l,m0 

Setup Sample buffer mod 

17 

move 

m0,m3 

Alternate Sample Buffer mod 

18 




19 

move 

M0, xO 


20 

do 

(POINTS , CLEAR.SAMPLE 

; Clear Sample Buffer 

21 

move 

x0,x:(r0) ;Real 

22 

move 

x0,y:(r0)+ ; Imaginary 

23 

CLEAR, SAMPLE 



24 




25 

move 

#FFT_DATA,rl 

Setup FFT Data ptr 

26 

move 

rl ,r7 

Imag. Pointer to FFT Buffer 

27 

move 

*POINTS-l,ml ;S«tup FFT Data mod 

28 

move 

ml ,m7 

Imag. FFT Buffer mod 

29 

do 

•POINTS. CLEAR.DATA 

Clear FFT Data buffer 

30 

move 

x0,x: (rl) 

Real 

31 

move 

x0,y : (rl)+ 

Imaginary 

32 

CLEAR.DATA 



33 




34 

; move 

•FFT.RESULT , r 2 

Setup FFT Result ptr 

35 

move 

•MAG.SQ.DATA , r2 


36 

move 

•POINTS/2-1 .m2 

Setup FFT Result mod 

37 

do 

•POINTS/2 . CLEAR.MAGSQ -.Clear FFT Result buffer 

38 

move 

xO.x: <r2) 

Real 

39 

move 

x0,y:(r2)+ 

Imaginary 

40 

CLEAR.MAGSQ 



41 




42 

move 

r7 ,x : IFFT.PTR 

Store IFFT for r7 reuse. 

43 

move 

m7 , x : IFFT.MOD 


44 

move 

#CNV_0UT,r7 

Setup Convolution Output ptr 

45 

move 

#2*P0INTS-2,m7 

Setup Conv. Output mod. 

46 

do 

#2»P0INTS-1 .CLEAR.CNV ;Clear Conv. Output 

47 

move 

xO.y: (r7)+ 


48 

CLEAR.CKV 



49 

move 

x: IFFT.PTR, r 7 


50 

move 

X : IFFT.MOD ,m7 


51 




52 

move 

*HAMM,r4 ; Setup Hamming ptr 

53 

move 

•POINTS-1 ,m4 

Setup Hamming mod 

54 

move 

#TAB,r5 

Sine Table 

55 

move 

tTABLE.SIZE-1 ,m5 

Sine Table mod 

56 




57 

move 

•POINTS, xl 

Initialize Sample Counter 

58 

move 

xl , x : SAMPLE.COUNTER 


59 

move 

•>ITERS,xl 

Init FFT Counter 

60 

move 

xl ,x:FFT.COUNTER 
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61 



62 

move 

xO,x: FLAGS 

63 



64 ; 

move 

•OUT.LENGTH , xl 

65 ; 

move 

xl ,x : OUT. COUNTER 

66 



67 

rts 



;Clear User Defined Flag Register 
;Set output duration counter. 


no 



hamming.asm 


01 hamming macro 

points 

tt 

i 

O 

n 

>— 

02 hamming idant 

1,2 



03 




04 py 

equ 

3.141592654 


06 FREQ. INC 

equ 

2.0»py/<cvf (pointa-1) frequency increment 

-■ 

06 SCALE.SHIFT 

equ 

6CVI(Clog(Ccvf (points))/Clog(2.0)) jShifts to Produce 1/P0INTS 


07 ; ITERATIONS 

equ 

2 ; Number of FFT Iterations to perform. There 


08 


; exists a limit before overflow. 


09 ; SCALE. FAC 

equ 

Ccvf (points)*#cvf (ITERATIONS) ;Scale Factor 


10 




11 org y : hamm. 

.loc 



12 N 

set 

0 


13 

dup 

points 


14 

dc 

(0. 54-0. 46»ecos(FREQ.INC*ecvf(N)))/6cvf (points/2) 


15 N 

set 

N+l 


16 

endm 



17 

endm 

;end of hamming macro 
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core-ipl.asm 


01 

02 

03 

04 

05 

06 

07 

08 

09 

10 
11 
12 

13 

14 
16 
16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 


****************************************************************************** 
Core Interrupt Priority Configuration Routine: 

Purpose: The purpose of this routine is to set the core Interrupt 

priorities. For the DAMA Project, only IRQB need concern us presently. 

Thus, only bits 5 to 3 are relevant. The following table suggests the 
proper settings: 

IBL2: 0 for level triggering, 1 for edge triggering (DAMA uses edge) 

IBL1-0: 

00 01 10 11 


Enabled 

Priority 


No Yes Yes Yes 
0 12 


For details see 303UM:D-17. 

Current Settings: 

Currently IRQB is the only interrupt enabled and it has been set to 
priority level 2 (highest) and negative edge triggering. 

Written By: Tim Bagget 

Adapted By: Brad Scaife 

Date: 3/22/98 

Notes: For use with rev 2.1 code. 

****************************************************************************** 
CORE Interrupt Priority and Configuration 
IBM equ tl ; IRQB trigger (0 level, 1 neg edge) 

IBP equ $1 ;IRQB priority level 0, 1, or 2; 

IBL equ (IBM«2) + (IBP+l)Jt3 

;IPRC equ IBL<<M_IBL0tM_IBL ; Disabled for the time being. 

IPRC equ $000038 


movep #IPRC,x:M_IPRC 


llnitialize Interrupt Priority/Conf ig 


112 



sintab.asm 


01 sintab 

02 sintab 

03 

04 pie 

05 TAB.INC 

06 


07 

08 N 

09 

10 

11 N 

12 

13 

14 


org 


macro tasiz,tab_loc 

ident 1,2 

equ 3.141592654 

equ 2.0*pie/Ccvf (taaiz) 

x:tab_loc 

set 0 

dup tasiz 

dc #sin(TAB_INC*#cvi (N))/2.0 

set N+l 

endm 

endm 

;end sine table generation macro 
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7819equ.asm 


1 BB.ADR 

2 BB7819.DR 


equ *0 ; DIP Switch Address SW1 ($0 - $3f) 

equ $FFFF80+BB_ADR ; ADS7819 Data Register 
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fftr2cn.asm 


001 

002 

003 

004 
006 
006 

007 

008 

009 

010 

011 

012 

013 

014 


061 


; This program originally available on the Motorola DSP bulletin board. 
; It is provided under a DISCLAMER OF WARRANTY available from 
; Motorola DSP Operation, 6601 Wm. Cannon Drive W., Austin, Tx., 78736. 

; Radix 2, In-Place, Decimation-In-Time FFT (fast). 

• 

; Last Update 18 Aug 88 Version 1.0 

t 

fftr2cn macro points, data, odata.coef 

fftr2cn ident 1,0 

Radix 2 Decimation in Time In-Place Fast Fourier Transform Routine 


016 

Complex input 

and output data 


016 

Real 

data 

in X memory 


017 

Imaginary data in Y memory 


018 

Normally ordered input data 


019 

Normally ordered output data 


020 

Coefficient 

lookup table 


021 

-Cosine values in X memory 


022 

-Sine 

values in Y memory 


023 





024 

Macro Call - 

fftr2cn points, data, odata, coef 

025 





026 

points 


number of points (16 

■32768, 

027 

data 


start of data buffer 


028 

odata 


start of output data 

buffer 

029 

coef 


start of sine/cosine 

table 

030 





031 

Alters Data 

ALU Registers 


032 

xl 

xO 

yl yO 


033 

a2 

al 

aO a 


034 

b2 

bl 

bO b 


035 





036 

Alters Address Registers 


037 

rO 

nO 

mO 


038 

rl 

nl 

ml 


039 


n2 



040 





041 

r4 

n4 

m4 


042 

r6 

n5 

mS 


043 

r7 

n7 

m7 


044 





045 

Alters Program Control Registers 


046 

PC 

ar 



047 





048 

Uses 6 locations 

on System Stack 


049 





050 

Latest Revision - 

18 Aug-88 



062 

move fdata.rO 

; initialize 

053 

move #points/4.n0 

; initialize 

064 

move n0,n4 

; 

055 

move n0,n7 

; initialize 

056 

move #points-l,m0 

; initialize 

057 

move mO.ml 

;for modulo 

058 

move m0,m4 


059 

move m0,m5 



offset 


060 
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061 ; Do first and second Radix 2 FFT passes, combined as 4-point butterflies 

062 ; 


063 

move 


x: (rO)+nO,xO 



064 

tfr 

xO ,a 

x: (rO)+nO,yl 



065 






066 

do 

nO, .twopass 



067 

tfr 

yl,b 

x: (rO)+nO,yO 



068 

add 

y0,a 

x: (rO) ,xl 


;ar+cr 

069 

add 

xl,b 

r0,r4 


;br+dr 

070 

add 

a,b 

(rO)+nO 


;ar ’=(ar+cr)+(br+dr) 

071 

subl 

b.a 

b.x: (rO)+nO 


;br ’=(ar+cr)-(br+dr) 

072 

tfr 

x0,a 

a.xO 

y: (rO) ,b 


073 

sub 

y0,a 


y : (r4)+n4,y0 

;ar-cr 

074 

sub 

yo.b 

xO,x : (rO) 


;bi-di 

075 

add 

a ,b 


y : (rO)+nO,xO 

; cr ’=(ar-cr)+(bi-di) 

076 

subl 

b, a 

b.x: (rO) 


;dr ’=(ar-cr)-(bi-di) 

077 

tfr 

xO,a 

a.xO 

y : (r4) ,b 


078 

add 

yO.a 


y : (rO)+nO,yO 

jbi+di 

079 

add 

yO,b 

xO.x: (rO)+nO 


;ai+ci 

080 

add 

b.a 


y: (r0)+, xO 

;ai’=(ai+ci)+(bi+di) 

081 

subl 

a,b 


a.y : (r4)+n4 

;bi ’=(ai+ci)-(bi+di) 

082 

tfr 

xO,a 


b.y : (r4)+n4 


083 

sub 

yO,a 

xl ,b 


;ai-ci 

084 

sub 

yi.b 

x: (rO)+nO,xO 


;dr-br 

085 

add 

a,b 

x: (rO)+nO,yl 


;ci'=(ai-ci)+(dr-br) 

086 

subl 

b.a 


b.y: (r4)+n4 

;di ’=(ai-ci)-(dr-br) 

087 

tfr 

xO,a 


a.y: (r4)+ 



088 .twopass 

089 ; 

090 ; Perform all next FFT passes except last pass with triple nested DO loop 

091 ; 


092 

move 

tpoints/8,nl ; initialize butterflies 

per group 


093 

move 

#4,n2 

.-initialize groups per 

pass 


094 

move 

#-l,m2 

; linear addressing for 

r2 


095 

move 

tO ,m7 

.-initialize C address modifier for 


096 



; reverse carry (bit-reversed) addressing 


097 






098 

do 

#Ocvi(®log(points)/01og(2)-2 .5) ..end.pass 

;-l ??? example: 7 passes for 1024 pt 

. FFT 

099 

move 

tdata.rO 


.•initialize A input pointer 


100 

move 

rO.rl 




101 

move 

nl ,r2 




102 

move 

rO ,r4 


; initialize A output pointer 


103 

move 

(rl)+nl 


; initialize 8 input pointer 


104 

move 

rl ,r5 


; initialize B output pointer 


105 

move 

tcoef ,r7 


; initialize C input pointer 


106 

lua 

(r2)+,n0 


initialize pointer offsets 


107 

move 

n0,n4 




108 

move 

n0,n5 




109 

move 

(r2)- 


; butterfly loop count 


110 

move 


x: (rl) ,xl y:(r7),y0 

;xl=br,yO=vi .lookup -sine and -cosine values 

111 

move 


x:(r7)+n7,x0 y:(rO),b 

; b=ai, xO=vr, update C pointer, preload data 

112 

mac 

xl.yO.b 

y:(rl)+,yl 

;yl=bi ,b=ai+br«wi. 


113 

macr 

-xO.yl.b 

y:(rO).a 

; a=ai again, b=ai+br*wi-bi *wr=bi ’ 


114 






115 

do 

n2 ..end.grp 



116 

do 

r2,_end.bfy 

;loop BF-1 times and start Radix 2 

DIT BF kernel 

117 

subl 

b.a 

x:(rO),b b,y:(r4) 

;b=ar ,a=ai-br*wi+bi*vr=ai ’ , 

PUT bi’ 

118 

mac 

-xl ,xO,b 

x:(rO)+,a a,y:(r5) 

; a=br . b=ar-br*wr , 

PUT ai’ 

119 

macr 

-yl ,yO,b 

x: (rl) ,xl 

;b*ar- (br»wr+bi*wi)=br ’ ,xl=next br 


120 

subl 

b.a 

b,x:(r4)+ y:(rO),b 

;b=nai ,a=2*ar-ar+br*wr+bi*wi=ar 1 . 

PUT br’ 

121 

mac 

xl.yO.b 

y : (rl)+,yl 

;yl»nbi ,b=nai+nbr*wi 


122 

macr 

-xO.yl ,b 

a,x:(r5)+ y:(rO),a 

;a=nai ,b=nai+nbr*wi-nbi»wr=nbi ’ 

PUT ar’ 
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123 .end.bfy 


124 

move 

(rl)+nl 



points to first B in next 

group 


125 

subl 

b.a 

x: (rO) ,b 

b.y: (r4) 

PUT last bi’ 

in a group 



126 

mac 

-xl.xO.b 

x: (r0)+n0,a 

a.y: (rS) 

PUT last ai’ 

in a group 

and update 

A pointer 

127 

macr 

-yl.yO.b 

x: (rl) ,xl 

y : (r7) ,y0 





128 

subl 

b.a 

b.x: (r4)+n4 

y: (rO) ,b 

PUT last br’ 

in a group 

and update 

A’ pointer 

129 

mac 

xl.yO.b 

x:(r7)+n7,x0 

y: (rl)+,yl 

update W pointer 



130 

macr 

-xO.yl.b 

a.x: (r5)+n5 

y: (rO) ,a 

PUT last ar' 

in a group 

and update 

B* pointer 

131 

_«nd_grp 








132 

move 

nl.bl 







133 

lsr 

b n2 , a: 

'.divide 

butterflies per group by two 





134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 


lsl 

move 

.end.pass 


bl,nl 


al , n2 


Do last FFT pass 

move #2,n0 
move nO.nl 
move *points/4,n4 
move n4,n5 
move tdata.rO 
move *odata,r4 
move r4,r2 
lua (rO)+,rl 
lua (r2)+n2,r5 
move #0,m4 
move m4,m5 
move #coef,r7 
move (r5)-n5 


; multiply groups per pass by two 


; initialize pointer offsets 

; output pointer A offset 
; output pointer B offset 
; initialize A input pointer 
; initialize A output pointer 
;save A output pointer 
; initialize B input pointer 
; initialize B output pointer 
;bit-reversed addressing for output ptr. A 
; bit-reversed addressing for output ptr. B 
; initialize C input pointer 
; predecrement output pointer 


153 

164 

155 

move 

move 


x: (rl) ,xl 
x: (r5) .a 

y:(r7),y0 ;xl*br,yO«wi 
y: (rO),b ;a«?,b>»ai 



166 

157 

do 

mac 

n2. .lastpass ; Radix 

xl.yO.b x:(r7)+n7,x0 

2 DIT butterfly kernel with one butterfly per group 
y:(rl)+nl,yl ;b«ai+br»wi,xO«wr , yl«bi 

158 

macr 

-xO.yl.b 

a.x: (r5)+n5 

y:(rO) ,a ;b«ai+br*wi-bi*wr=bi’ ,a*ai. 

PUT 

previous ar 

159 

subl 

b ,a 

x: (rO) ,b 

b,y:(r4) ;a=ai’,b»ar. 

PUT 

bi* 

160 

161 

mac 

macr 

-xl.xO.b 

-yl.yO.b 

x : (rO)+nO,a 
x: (rl) ,xl 

a,y:(r5) ; b«ar-br»wr ,a«ar , 

y:(r7),y0 ;b»br ’ ,xl=nbr ,yO=nwi 

PUT 

ai’ 

162 

163 

subl 

.lastpass 

b.a 

b.x: (r4)+n4 

y:(rO),b ;a*ar’ ,b=nai. 

PUT 

br ’ 

164 

165 

move 

endm 


a.x: (r5)+n5 

: 

PUT 

ar' 
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sincos.asm 


01 

02 

03 

04 

05 

06 

07 

08 

09 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 


; This program originally available on the Motorola DSP bulletin board 
; It is provided under a DISCLAMER OF WARRANTY available from 
; Motorola DSP Operation, 6501 Wm. Cannon Drive W. , Austin, Tx., 78735 

; Sine-Cosine Table Generator for FFTs. 

; Last Update 25 Nov 86 Version 1.2 

sincos macro points, coef 

sincos ident 1,2 

sincos - macro to generate sine and cosine coefficient 

lookup tables for Decimation in Time FFT 
twiddle factors. 

points - number of points (2 - 32768, power of 2) 

coef - base address of sine/cosine table 

negative cosine value in X memory 
negative sine value in Y memory 

Latest revision - 25-Nov-86 


25 pi 

equ 

3.141592654 

26 freq 

equ 

2. 0*pi/Ccvf (points) 

27 



28 

org 

x:coef 

29 count 

set 

0 

30 

dup 

points/2 

31 

dc 

-6cos(Ccvf (count )*freq) 

32 count 

set 

count+1 

33 

endm 


34 



35 

org 

y :coef 

36 count 

set 

0 

37 

dup 

points/2 

38 

dc 

-#sin(#cvf (count )*f req) 

39 count 

Bet 

count+1 

40 

endm 


41 



42 

endm 

;end of sincos macro 
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ada_equ.asm 


001 

002 

003 

004 

005 

006 

007 

008 

009 

010 

011 

012 

013 

014 

015 

016 

017 

018 

019 

020 
021 
022 

023 

024 


page 132,60 

********************************* ************************* ****************** 
ADA.EQU. ASM 

Initialization constants to facilitate initialization of the CS4215 

Copyright (c) MOTOROLA 1996 

Semiconductor Products Sector 
Digital Signal Processing Division 

**************************************************************************** 


025 N0.PREAMP 

026 

equ 

$100000 



027 LO.OUT.DRV 

028 

equ 

$080000 



029 HI_PASS_FILT 

030 

equ 

$008000 



031 SAMP.RATE.9 

032 

equ 

$003800 

; 9.6 kHz sample rata 


033 SAMP_RATE_48 

equ 

$003000 

; 48 kHz sample rate 


034 

036 SAMP_RATE_32 
036 

equ 

$001800 

; 32 kHz sample rate 

— 

037 SAMP.RATE.27 

038 

equ 

$001000 



039 SAHP.RATE.16 

040 

equ 

$000800 


— 

041 SAMP_RATE_8 

042 

equ 

$000000 



043 STEREO 

equ 

$000400 



044 

045 DATA.8LIN 

046 

equ 

$200300 



047 DATA.8A 

equ 

$200200 



048 

049 DATA.8U 

050 

equ 

$200100 



051 DATA.16 

052 

equ 

$200000 



053 IHMED.3STATE 

equ 

$800000 



054 

055 XTALl.SELECT 

056 

equ 

$100000 

; 24.576 MHz 


057 XTAL2. SELECT 

058 

equ 

$200000 

; 16.9344 MHz 


059 BITS.64 

060 

equ 

$000000 
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$040000 


061 BITS.128 

062 

063 BITS.256 

064 

065 CODEC.MASTER 

066 

067 CODEC.TX.OFF 

068 

069 

070 

071 ;CTRl.WD_12 

072 

073 ;CTRL_WD.34 

074 

075 ;CTRL_WD_56 

076 

077 ;CTRL.WD_78 

078 

079 

080 

081 HEADPHONE.EN 

082 

083 LINED UT.EN 

084 

085 SPEAKER.EN 

086 

087 MIC.IN.SELECT 

088 

089 LEFT.ATTN 

090 

091 RIGHT.ATTN 

092 

093 LEFT.GAIN 

094 

095 RIGHT.GAIN 

096 

097 MONITOR_ATTN 

098 

099 ; OUTPUT_SET 

100 

101 ; INPUT.SET 


equ 

equ $080000 

equ $020000 

equ $010000 

equ N0_PREAMP+HI_PASS_FILT+SAMP_RATE_48+STERE0+DATA_16 
equ IMMED_3STATE+XTAL1_SELECT+BITS_64+C0DEC. MASTER 
equ $000000 

equ $000000 

equ $800000 

equ $400000 

equ $004000 

equ $100000 


equ 

$010000 

;63*LEFT_ATTN 

-94.5 

dB, 

1.5 

dB 

steps 

equ 

$000100 

;63*RIGHT.ATTN = 

-94.5 

dB, 

1.5 

dB 

steps 

equ 

$010000 

; 15»LEFT_GAIN 

22.5 

dB, 

1.5 

dB 

steps 

equ 

$000100 

;16*RIGHT.GAIN 

22.5 

dB, 

1.5 

dB 

steps 

equ 

$001000 

; 15*M0NIT0R_ATTN = 

mute , 


6 

dB 

steps 


equ HEADPH0NE_EN+LINE0UT_EN+(LEFT_ATTN»4) 

equ MIC. IN.SELECT + ( 15*MQNIT0R_ ATTN) + (RIGHT. ATTN *4) 


;CLB=0 
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vectors, asm 


001 

002 

003 

004 

005 

006 

007 

008 

009 

010 
Oil 


page 132,60 

**************************************************************************** 
VECTORS. ASM 

Vector table for the 56303 

Copyright (c) MOTOROLA 1996 

Semiconductor Products Sector 
Digital Signal Processing Division 

*************** ******************************** ************ ******** •**•**••• 


; Hardware RESET 


012 

013 ORG P:0 

014 ; 

015 vectors JMP START 

016 ; 

017 jmp * 

018 NOP 

019 

020 jmp * 

021 NOP 

022 ; 

023 jmp * 

024 NOP 

025 ; 

026 jmp * 

027 NOP 

028 ; 

029 jmp * 

030 NOP 

031 ; 

032 NOP 

033 NOP 

034 ; 

035 NOP 

036 NOP 

037 ; 

038 jsr main 

039 ; 

040 jmp * 

041 NOP 

042 ; 

043 jmp * 

044 NOP 

045 ; 

046 jsr echo 

047 ; 

048 jmp * 

049 NOP 

050 ; 

061 jmp * 

052 NOP 

053 ; 

054 jmp * 

055 NOP 

056 ; 

057 jmp * 

058 NOP 

059 ; 

060 jmp * 


Stack Error 
Debug Request Interrupt 
Debug Request Interrupt 
Trap 
NMI 

Reserved 

Reserved 

IRQA 

IRQB 

IRQC 

IRqD 

DMA Channel 0 
DMA Channel 1 
DMA Channel 2 
DMA Channel 3 
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061 

NOP 


; - DMA Channel 4 

062 ; 




063 

jmp 

* 


064 

NOP 


DMA Channel 5 

065 ; 




066 

jmp 

* 


067 

NOP 


Timer 0 Compare 

068 ; 




069 

j“P 

* 


070 

NOP 


Timer 0 Overflow 

071 ; 




072 

jmp 

* 


073 

NOP 


Timer 1 Compare 

074 ; 




075 

j“P 

* 


076 

NOP 


Timer 1 Overflow 

077 ; 




078 

jmp 

* 


079 

NOP 


Timer 2 Compare 

080 ; 




081 

jmp 

* 


082 

NOP 


Timer 2 Overflow 

083 ; 




084 

jsr 

ssi_rx_isr 

ESSIO Receive Data 

085 ; 




086 

jsr 

ssi.rxe.isr 

ESSIO Receive Data w/ Exception Status 

087 ; 




088 

jsr 

ssi_rxls_isr 

ESSIO Receive last slot 

089 ; 




090 

jsr 

ssi.tx.isr 

ESSIO Transmit Data 

091 ; 




092 

jsr 

ssi._txe.isr 

ESSIO Transmit Data w / Exception Status 

093 ; 




094 

jsr 

ssi_txls_isr 

ESSIO Transmit last slot 

095 ; 




096 

NOP 


Reserved 

097 

NOP 



098 ; 




099 

NOP 


Reserved 

100 

NOP 



101 ; 




102 

jmp 

* 


103 

NOP 


ESSI1 Receive Data 

104 ; 




105 

jmp 

* 


106 

NOP 


ESSI1 Receive Data w/ Exception Status 

107 ; 




108 

jmp 

* 


109 

NOP 


ESSI1 Receive last slot 

110 ; 




111 

jmp 

* 


112 

NOP 


ESSI1 Transmit Data 

113 ; 




114 

jmp 

* 


115 

NOP 


ESSI1 Transmit Data w/ Exception Status 

116 ; 




117 

jmp 

* 


118 

NOP 


ESSI1 Transmit last slot 

119 ; 




120 




121 

NOP 


Reserved 

122 

NOP 




122 



185 

jmp * 





186 

NOP 

; Available 

for 

Host 

Command 

187 

j B p * 





188 

NOP 

; Available 

for 

Host 

Command 

189 

jnp * 





190 

NOP 

; Available 

for 

Host 

Command 

191 

jmp * 





192 

NOP 

; Available 

for 

Host 

Command 

193 

jop * 





194 

NOP 

; Available 

for 

Host 

Command 

195 

jmp « 





196 

NOP 

; Available 

for 

Host 

Command 

197 

jmp * 





198 

NOP 

; Available 

for 

Host 

Command 

199 

jmp » 





200 

NOP 

; Available 

for 

Host 

Command 

201 

jmp * 





202 

NOP 

; Available 

for 

Host 

Command 

203 

jmp * 





204 

NOP 

; Available 

for 

Host 

Command 

205 

jmp » 





206 

NOP 

; Available 

for 

Host 

Command 

207 

jmp » 





208 

NOP 

; Available 

for 

Host 

Command 

209 

jmp * 





210 

NOP 

; Available 

for 

Host 

Command 

211 

jmp * 





212 

NOP 

; Available 

for 

Host 

Command 

213 

jmp » 





214 

NOP 

; Available 

for 

Host 

Command 

215 

jmp * 





216 

NOP 

; Available 

for 

Host 

Command 

217 

jmp « 





218 

NOP 

; Available 

for 

Host 

Command 

219 

jmp * 





220 






221 

NOP 

; Available 

for 

Host 

Command 

222 

jmp » 





223 

NOP 

; Available 

for 

Host 

Command 

224 

jmp » 





225 

NOP 

; Available 

for 

Host 

Command 

226 

jmp * 





227 

NOP 

I Available 

for 

Host 

Command 

228 

jmp « 





229 

NOP 

; Available 

for 

Host 

Command 

230 

jmp » 





231 

NOP 

; Available 

for 

Host 

Command 

232 

jmp * 





233 

NOP 

; Available 

for 

Host 

Command 

234 

jmp * 





236 

NOP 

; Available 

for 

Host 

Command 

236 

jmp * 





237 

NOP 

; Available 

for 

Host 

Command 

238 

jmp • 





239 

NOP 

; Available 

for 

Host 

Command 

240 

jmp * 





241 

NOP 

; Available 

for 

Host 

Command 

242 

jmp ♦ 





243 

NOP 

; Available 

for 

Host 

Command 

244 

jmp * 





245 

NOP 

; Available 

for 

Host 

Command 

246 

jmp » 
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247 

NOP 


; Available 

for 

Host 

Command 

248 

jmp 

* 





249 

NOP 


; Available 

for 

Host 

Command 

260 

jmp 

4 





251 

NOP 


; Available 

for 

Host 

Command 

252 

jmp 

4 





263 

NOP 


; Available 

for 

Host 

Command 

254 

jmp 

4 





255 

NOP 


; Available 

for 

Host 

Command 

256 

jmp 

4 





257 

NOP 


; Available 

for 

Host 

Command 

258 

jmp 

4 





259 

NOP 


; Available 

for 

Host 

Command 

260 

jmp 

4 





261 

NOP 


; Available 

for 

Host 

Command 

262 

jmp 

4 





263 

NOP 


; Available 

for 

Host 

Command 

264 

jmp 

4 





265 

NOP 


; Available 

for 

Host 

Command 

266 

jap 

4 





267 

NOP 


; Available 

for 

Host 

Command 

268 

jmp 

4 





269 

NOP 


; Available 

for 

Host 

Command 

270 

jmp 

4 





271 

NOP 


; Available 

for 

Host 

Command 

272 

jmp 

4 





273 

NOP 


; Available 

for 

Host 

Command 

274 

jmp 

4 





275 

NOP 


; Available 

for 

Host 

Command 

276 

jmp 

4 





277 

NOP 


; Available 

for 

Host 

Command 

278 

jmp 

4 





279 

NOP 


; Available 

for 

Host 

Command 

280 

jap 

4 





281 

NOP 


; Available 

for 

Host 

Command 

282 

jmp 

4 





283 







284 







285 

NOP 


; Available 

for 

Host 

Command 

286 

jmp 

4 





287 

NOP 


; Available 

for 

Host 

Command 

288 

jmp 

4 





289 

NOP 


; Available 

for 

Host 

Command 

290 

jmp 

4 





291 

NOP 


; Available 

for 

Host 

Command 

292 

jmp 

4 





293 

NOP 


; Available 

for 

Host 

Command 

294 

jap 

4 





295 

NOP 


; Available 

for 

Host 

Command 

296 

jmp 

4 





297 

NOP 


; Available 

for 

Host 

Command 

298 

jmp 

4 





299 

NOP 


; Available 

for 

Host 

Command 

300 

jmp 

4 





301 

NOP 


; Available 

for 

Host 

Command 

302 

jmp 

4 





303 

NOP 


; Available 

for 

Host 

Command 

304 

jmp 

4 





305 

NOP 


; Available 

for 

Host 

Command 

306 

jap 

4 





307 

NOP 


; Available 

for 

Host 

Command 

308 

jap 

4 
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309 

NOP 


310 

jmp 

* 

311 

NOP 


312 

jmp 

* 

313 

NOP 


314 

jmp 

* 

315 

NOP 


316 

jmp 

* 

317 

NOP 


318 ; 



319 ; 




Available 

for 

Host 

Command 

Available 

for 

Host 

Command 

Available 

for 

Host 

Command 

Available 

for 

Host 

Command 

Available 

for 

Host 

Command 
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intequ.asm 


01 

02 

03 

04 
06 
06 

07 

08 

09 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 
63 

54 

55 

56 

57 

58 

59 

60 


****************** 4 *************** ******************************* **************** 

EQUATES for ONYXE 56302 interrupts 
Last update: June 11 1995 

44 * 4 ***** 44 ************************ 44 ******************************************** 

page 132,55,0,0,0 
opt mex 

intequ ident 1,0 

if CDEFCI.VEC) 

; leave user definition as is. 
else 

I.VEC equ $0 
endif 

; Non-Maskable interrupts 


.RESET 

EQU 

I_VEC+*00 

Hardware RESET 

.STACK 

Equ 

I.VEC+J02 

Stack Error 

.ILL 

EQU 

I.VEC+J04 

Illegal Instruction 

.DBG 

EQU 

I.VEC+S06 

Debug Request 

.TRAP 

EQU 

I.VEC+108 

Trap 

_NMI 

EQU 

I_VEC+t0A 

Non Maskable Interrupt 

Interrupt Request Pins 

-IRQA 

EQU 

I.VEC+S10 

IRQA 

_IRQB 

EQU 

I.VEC+112 

IRQB 

_IRQC 

EQU 

I.VEC+114 

IRQC 

.IRQD 

EQU 

I.VEC+116 

IRQD 


DMA Interrupts 


t.DMAO 

EQU 

I.VEC+118 

DMA Channel 0 

[.DMA1 

EQU 

I_VEC+*1A 

DMA Channel 1 

[_DMA2 

EQU 

I.VEC+$1C 

DMA Channel 2 

[.DM A3 

EQU 

I.VEC+J1E 

DMA Channel 3 

I.DMA4 

EQU 

I.VEC+J20 

DMA Channel 4 

_DMA5 

EQU 

I.VEC+*22 

DMA Channel 5 

Timer 

Interrupts 


[.TIMOC 

EQU 

I.VEC+124 

TIMER 0 compare 

[.TIMOOF EQU 

I.VEC+$26 

TIMER 0 overflow 

[.TIM1C 

EQU 

I.VEC+»28 

TIMER 1 compare 

I.TIM10F 

EQU 

I.VEC+J2A 

TIMER 1 overflow 

I.TIM2C 

EQU 

I.VEC+I2C 

TIMER 2 compare 

[.TIM20F 

EQU 

I.VEC+I2E 

TIMER 2 overflow 

ESSI Interrupts 

[.SIORD 

EQU 

I.VEC+130 

ESSIO Receive Data 
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61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 


I.SIORDE 

EQU 

I_VEC+$32 

ESSIO Receive Data With Exception 

Status 

I.SIORLS 

EQU 

I.VEC+634 

ESSIO Receive last slot 


I.SIOTD 

EQU 

I.VEC+136 

ESSIO Transmit data 


I.SIOTDE 

EQU 

I.VEC+138 

ESSIO Transmit Data With Exception 

Status 

I.SIOTLS 

EQU 

I.VEC+I3A 

ESSIO Transmit last slot 


I. SURD 

EQU 

I.VEC+140 

ESSI1 Receive Data 


I.SI1RDE 

EQU 

I_VEC+$42 

ESSI1 Receive Data With Exception 

Status 

I.SIIRLS 

EQU 

I_VEC+*44 

ESSI1 Receive last slot 


I.SI1TD 

EQU 

I.VEC+146 

ESSI1 Transmit data 


I.SI1TDE 

EQU 

I_VEC+$48 

ESSI1 Transmit Data With Exception 

Status 

I.SI1TLS 

EQU 

I_VEC+*4A 

ESSI1 Transmit last slot 


; SCI Interrupts 

I.SCIRD 

EQU 

I_VEC+*50 

SCI Receive Data 


I.SCIRDE 

EQU 

I.VEC+»52 

SCI Receive Data With Exception Status 

I.SCITD 

EQU 

I.VEC+154 

SCI Transmit Data 


I.SCIIL 

EQU 

I.VEC+166 

SCI Idle Line 


I.SCITM 

EQU 

I.VEC+658 

SCI Timer 


; HOST Interrupts 

I.HRDF 

EQU 

I.VEC+160 

; Host Receive Data Full 


I.HTDE 

EQU 

I.VEC+162 

; Host Transmit Data Empty 


I.HC 

EQU 

I.VEC+164 

; Default Host Command 


; INTERRUPT ENDING ADDRESS 

I. INTEND 

EQU 

I_VEC+$FF 

last address of interrupt vector space 


LIST 
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sincos.asm 


01 ; 

02 ; This program originally available on the Motorola DSP bulletin board 

03 ; It is provided under a DISCLAMER OF WARRANTY available from 

04 ; Motorola DSP Operation, 6501 Wm. Cannon Drive W., Austin, Tx., 78735 

05 ; 

06 ; Sine-Cosine Table Generator for FFTs . 

07 ; 

08 ; Last Update 25 Nov 86 Version 1.2 


09 




10 sincos 

macro 

points , coef 

11 sincos 

ident 

1.2 

12 




13 


sincos 

macro to generate sine and cosine coeffi 

14 



lookup tables for Decimation in Time FFT 

15 



twiddle factors. 

16 




17 


points 

number of points (2 - 32768, power of 2) 

18 


coef 

base address of sine/cosine table 

19 



negative cosine value in X memory 

20 



negative sine value in Y memory 

21 




22 

Latest revision - 25-Nov-86 

23 




24 




25 pi 

equ 

3.141592654 

26 freq 

equ 

2. 0*pi/«cvf (points) 

27 




28 


org 

x:coef 

29 count 

set 

0 

30 


dup 

points/2 

31 


dc 

-®cos(®cvf (count )*freq) 

32 count 

set 

count+1 

33 


endm 


34 




35 


org 

y : coef 

36 count 

set 

0 

37 


dup 

points/2 

38 


dc 

-®sin(8cvf (count )*f req) 

39 count 

set 

count+1 

40 


endm 


41 




42 


endm 

;end of sincos macro 
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adaJnit.asm 


001 

002 

003 

004 

005 

006 

007 

008 

009 

010 

011 

012 

013 

014 

015 

016 

017 

018 

019 

020 
021 
022 

023 

024 

025 

026 

027 

028 

029 

030 

031 

032 

033 

034 

035 

036 

037 

038 

039 

040 

041 

042 

043 

044 

045 

046 


page 132,60 

**************************************************************************** 
ADA.INIT. ASM Var.2.0 
Example program to initialize the CS4215 

Copyright (c) MOTOROLA 1995, 1996 

Semiconductor Products Sector 
Digital Signal Processing Division 

History: 

14 June 1996: RLR/LJD - ver.1.0 

*************************************************************************** 

PLEASE NOTE: For use with rev 2.1 code. 
************************************************************************** 

portc usage: 

bit8: SSI TX (from DSP to Codec) 

bit7 : 

bit6: 

bitB: 

bit4: codec reset (from DSP to Codec) 
bit3: 

bit2: data/control bar 
0=control 
1-data 

«•******••* •••*•*•******** *********** ********** ********* ********* ********** 
***** initialize the CS4215 codec ***** 

****•••*•*••**•••••***•*•*•*•**••••*•••••**•••**•*****•••*******•**••**•**• 


PROGRAM OUTLINE: 

1 program fsync and sclk ** output 

2 write pcO * 0 (control mode) 

3 send 64 bit frame x times, with deb bit * 0, keep doing until read back as 0 

4 send 64 bit frame x times, with deb bit = 1, keep doing until read back as 1 

5 re-program fsync and sclk »» input 

6 write pcO « 1 (data mode) 

7 receive/send data (echo slots 1,2, 3, 4; slots 5, 6, 7, 8 == constants) 
**************************************************************************** 


initialize ssi — fsync and sclk ==> outputs 


047 

org 

p ; 


048 

ada_init 



049 

movep 

MOOOO.x :M_PCRC 

; turn off ESSIO port (for now) 

050 

movep 

#1103807 ,x:M_CRA0 

; 40MHz/16 = 2.5MHz SCLK, WL=16 bits, 4H/F 

051 

movep 

#tf f 313C , x : M_CRB0 

; RIE,TIE,RLIE,TLIE,RE,TE, sc2/sck outputs 

062 

movep 

#10003 , x : M_PRRC 

; setup pdO and pdl as gpio output 

053 

movep 

#*0,x:M.PDRC 

; send out a 0 on DC" and RST.CODEC" 

054 




055 

; — 

— reset delay for codec 

— 

056 

do 

#1000,_delay_loop 


057 

rep 

#2000 

; 100 us delay (assuming 40MHz VCO) 

058 

nop 



059 

_delay_loop 



060 
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061 

bset 

#0,x:M.PDRC 

; sends out a 1 on pdO (rst_codec*l) 

062 

movep 

#10008, x;M.IPRP 

; set interrupt priority level for ESSIO to 

063 

andi 

#*FC,mr 


; enable interrupts 

064 





066 

***************************************************************************** 

066 

The following 

data sets up the CS4215 

control mode data: 

067 

(CTS = Control Time Slot, U/LN = upper/lower Nibble) 

068 





069 

+ 

CTS 1 -UN 

0 0 

1 MLB 0000 

070 

1+ 

CTS1-LN 

□LB CLB 

X X 0 0 0 0 

071 

11+ 

CTS2-UN 

HPF X 

DFR2 DFR1 0 0 10 

072 

111+ 

CTS2-LN 

DFRO ST 

DF1 DFO 1100 

073 

xO = $002Cxx 




074 





076 

+ 

CTS3-UN 

ITS MCK2 

MCK1 MCKO 1000 

076 

|+ 

CTS3-LN 

BSEL1 BSELO 

XCLK XEN 1000 

077 

II+- 

CTS4-UN 

TEST TEST 

TEST TEST (TEST MUST BE 0) 

078 

111+ 

CTS4-LN 

TEST TEST 

ENL DAD 0000 

079 

xO ■ $8800xx 




080 

***************************************************************************** 

081 





082 

set up buffer with 

control mode data 

083 

move 

#CTRL_WD.12,xO 


084 

move 

xO,x:TX 

.BUFF.BASE 


085 

move 

#CTRL_WD_34 ,x0 


086 

move 

xO.x.TX. 

BUFF.BASE+l 


087 

move 

#CTRL.WD_56,xO 


088 

move 

xO , x : TX. 

.BUFF.BASE+2 


089 

move 

#CTRL_WD.78,xO 


090 

move 

xO,x:TX. 

BUFF.BASE+3 


091 





092 

movep 

#*003C,x:M_PCRC 

[turn on ESSIO except for scO and sc2 

093 





094 





095 

CLB »« 0 




096 





097 

jclr 

#3,x:M_SSISR0,* 

; wait until rx frame bit==l 

098 

jset 

#3,x :M_SSISR0 ,* 

; wait until rx frame bit==0 

099 

jclr 

#3 , x : M.SSISRO , » 

; wait until rx frame bit==l 

100 

jset 

#18 ,x:RX_BUFF_BASE , * 

; loop until CLB set 

101 





102 





103 

CLB «= 1 




104 





105 

bset 

#18 , x : TX.BUFF.BASE 

;set CLB 

106 

do 

#4, .init.loopB 


107 

jclr 

#2 ,x :M_SSISR0 , * 

; wait until tx frame bit==l 

108 

jset 

#2,x:M_SSISR0,» 

; wait until tx frame bit==0 

109 

init.loopB 




110 

movep 

#30000, x:M_PCRC 

; disable ESSIO 

111 





112 

******* ******************* ****************************************** ********* 

113 

now CLB should be 1 — re-program fsync and sclk direction (i/p) — also, 

114 

circular buffer pointers for echoing data rO*current, rl=old data to send 

115 

1 frame later 



116 





117 

movep 

#3 103807, x:M_CRA0 

; 40MHz/16 * 2.5MHz SCLK, WL=16 bits, 4W/F 

118 

movep 

#3FF31OC,x:M.CRB0 

; sckd and fsync (sc02) as inputs 

119 

movep 

#30003, j 

. :M_PDRC 

; D/C" pin = 1 ==> data mode 

120 

movep 

#$003C , x : M.PCRC 

; turn on ESSIO except for scO and sc2 

121 

rts 




122 
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123 

**************************************************************************** 

124 

SSIO.ISR.ASM Ver .2.0 


125 

Example program to handle interrupts through 

126 

the 56303 SSI0 to move audio through the CS4215 

127 




128 

Copyright 

(c) MOTOROLA 1995, 1996 


129 


Semiconductor Products Sector 

130 


Digital Signal Processing 

Division 

131 




132 

upon 

entry : 


133 


R6 must be the stack pointer 

134 

corrupts: 


135 


R6 


136 




137 

History: 


138 


14 June 1996: RLR/LJD - 

ver 1.0 

139 

***************************************************************** ************* 

140 




141 




142 

the actual interrupt service routines (ISRs) follow: 

143 




144 

***•**«*•*•***•*****«•*• SSI TRANSMIT 

ISR ******************** • ******<■***»» 

145 ssi.txe.isr 



146 

bclr 

#4 , x : M.SSISR0 

; Read SSISR to clear exception flag 

147 



; explicitly clears underrun flag 

148 ssi.tx.isr 



149 

move 

r0,x: (r6)+ 

; Save rO to the stack. 

ISO 

move 

m0,x: (r6)+ 

; Save mO to the stack. 

151 

move 

*3 ,m0 

; Modulus 4 buffer. 

152 

move 

x:TX_PTR,r0 

; Load the pointer to the tx buffer. 

153 

nop 



154 

movep 

x: (rO)+,x:M_TX00 

; SSI transfer data register. 

155 

move 

rO , x : TX.PTR 

; Update tx buffer pointer. 

156 

move 

x:-(r6) ,m0 

; Restore mO. 

157 

move 

x:-(r6),r0 

; Restore rO. 

158 

rti 



159 




160 

********************* SSI TRANSMIT LAST SLOT ISR ************************** 

161 ssi.txls.isr 



162 

move 

r0,x: (r6)+ 

; Save rO to the stack. 

163 

move 

#TX_BUFF_BASE,r0 

; Reset pointer. 

164 

move 

r0,x:TX_PTR 

; Reset tx buffer pointer just in 

165 



; case it was corrupted . 

166 

move 

x:-(r6) ,r0 

; Restore rO. 

167 

rti 



168 




169 

************************** SSI receive 

ISR **••***•*•*•* ************ ******* 

170 ssi.rxe.isr 



171 

bclr 

#5.x:M_SSISR0 

; Read SSISR to clear exception flag 

172 



; explicitly clears overrun flag 

173 ssi.rx.isr 



174 

move 

r0,x : (r6)+ 

; Save rO to the stack. 

176 

move 

rnO.x: (r6)+ 

; Save mO to the stack. 

176 

move 

*3 ,m0 

; Modulo 4 buffer. 

177 

move 

x :RX_PTR,r0 

; Load the pointer to the rx buffer. 

178 

nop 



179 

movep 

x:M_RX0,x: (r0)+ 

; Read out received data to buffer. 

180 

move 

r0,x:RX.PTR 

; Update rx buffer pointer. 

181 

move 

x:-(r6) ,m0 

; Restore mO. 

182 

move 

x:-(r6) ,r0 

; Restore rO. 

183 

rti 



184 
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185 : ***************'******* SSI receive last slot ISR ************************** 


186 

ssi_rxls_isr 


187 

move 

r0,x: (r6)+ 

188 

move 

#RI. BUFF. BASE , rO 

189 



190 

move 

rO,x:RX.PTR 

191 

move 

x:-(r6) ,r0 

192 

rti 



; Save rO to the stack. 

; Reset rx buffer pointer just in 
; case it was corrupted. 

; Update rx buffer pointer. 

; Restore rO. 
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